The QCD Equation of State to \mathcal{O}(μ_B^6) from Lattice QCD
A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa, S. Mukherjee, H. Ohno, P. Petreczky, H. Sandmeyer, P. Steinbrecher, C. Schmidt, S. Sharma, W. Soeldner, M. Wagner
TThe QCD Equation of State to O ( µ B ) from Lattice QCD A. Bazavov, H.-T. Ding, P. Hegde a , O. Kaczmarek,
2, 4
F. Karsch,
4, 5
E. Laermann, Y. Maezawa, Swagato Mukherjee, H. Ohno,
5, 7
P. Petreczky, H. Sandmeyer, P. Steinbrecher,
4, 5
C. Schmidt, S. Sharma, W. Soeldner, and M. Wagner Department of Computational Mathematics, Science and Engineering and Department of Physics and Astronomy,Michigan State University, East Lansing, MI 48824, USA Key Laboratory of Quark & Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China Center for High Energy Physics, Indian Institute of Science, Bangalore 560012, India Fakult¨at f¨ur Physik, Universit¨at Bielefeld, D-33615 Bielefeld, Germany Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8317, Japan Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan Institut f¨ur Theoretische Physik, Universit¨at Regensburg, D-93040 Regensburg, Germany NVIDIA GmbH, D-52146 W¨urselen, Germany
We calculated the QCD equation of state using Taylor expansions that include contributions fromup to sixth order in the baryon, strangeness and electric charge chemical potentials. Calculationshave been performed with the Highly Improved Staggered Quark action in the temperature range T ∈ [135 MeV ,
330 MeV] using up to four different sets of lattice cut-offs corresponding to latticesof size N σ × N τ with aspect ratio N σ /N τ = 4 and N τ = 6 −
16. The strange quark mass is tunedto its physical value and we use two strange to light quark mass ratios m s /m l = 20 and 27, whichin the continuum limit correspond to a pion mass of about 160 MeV and 140 MeV respectively.Sixth-order results for Taylor expansion coefficients are used to estimate truncation errors of thefourth-order expansion. We show that truncation errors are small for baryon chemical potentialsless then twice the temperature ( µ B ≤ T ). The fourth-order equation of state thus is suitablefor the modeling of dense matter created in heavy ion collisions with center-of-mass energies downto √ s NN ∼
12 GeV. We provide a parametrization of basic thermodynamic quantities that canbe readily used in hydrodynamic simulation codes. The results on up to sixth order expansioncoefficients of bulk thermodynamics are used for the calculation of lines of constant pressure, energyand entropy densities in the T - µ B plane and are compared with the crossover line for the QCD chiraltransition as well as with experimental results on freeze-out parameters in heavy ion collisions. Thesecoefficients also provide estimates for the location of a possible critical point. We argue that resultson sixth order expansion coefficients disfavor the existence of a critical point in the QCD phasediagram for µ B /T ≤ T /T c ( µ B = 0) > . March 10, 2017
PACS numbers: 11.10.Wx, 12.38.Gc, 12.38Mh a [email protected] a r X i v : . [ h e p - l a t ] M a r I. INTRODUCTION
The temperature and density dependence of bulk thermodynamic quantities, commonly summarized as the equationof state (EoS), provide the most basic characterization of equilibrium properties of strong-interaction matter. Itsanalysis within the framework of lattice regularized Quantum Chromodynamics (QCD) has been refined ever sincethe early calculations performed in pure SU ( N ) gauge theories [1]. Quite recently, the continuum extrapolated resultsfor the EoS of QCD with physical light and strange quark masses have been calculated [2, 3]. Bulk thermodynamicobservables such as pressure ( P ), energy density ( (cid:15) ) and entropy density ( s ) as well as second order quantities such asthe specific heat ( C V ) and velocity of sound ( c s ) have now been obtained at vanishing chemical potentials for the threequark flavors ( µ u , µ d , µ s ). In accordance with the analysis of the chiral transition temperature, T c (cid:39) (154 ±
9) MeV [4],bulk thermodynamic observables change smoothly in the transition region. At low temperature they are found to bein quite good agreement with hadron resonance gas (HRG) model calculations, although some systematic deviationshave been observed, which may be attributed to the existence of additional resonances which are not taken intoaccount in HRG model calculations based on well established resonances listed in the particle data tables [5, 6].The EoS at vanishing chemical potentials does already provide important input into the modelling of the hy-drodynamic evolution of hot and dense matter created in heavy ion collisions. While this is appropriate for thethermal conditions met in these collisions at the LHC and the highest RHIC beam energies, knowledge of the EoSat non-vanishing baryon ( µ B ), strangeness ( µ S ) and electric charge ( µ Q ) chemical potentials is indispensable for thehydrodynamic modelling of the conditions met in the beam energy scan (BES) at RHIC. Due to the well-knownsign problem for lattice QCD formulations at non-zero chemical potential a direct calculation of the EoS at non-zero( µ B , µ Q , µ S ) is unfortunately not yet possible. At least for small values of the chemical potentials this can becircumvented by using a Taylor expansion of the thermodynamic potential [7, 8]. In this way some results for EoS atnon-zero baryon chemical potential have been obtained on coarse lattices [8–10]. These calculations have even beenextended to sixth order in the baryon chemical potential [11, 12]. First continuum extrapolated results for the EoSusing second order Taylor expansion coefficients have been obtained within the stout discretization scheme for stag-gered fermions [13] and simulations at imaginary chemical potential have been used to arrive at a sixth order resultfor the QCD EoS [14] and up to eighth order for some generalized susceptibilities [15] through analytic continuation.Results for higher order expansion coefficients are clearly needed if one wants to cover the range of chemicalpotentials, 0 ≤ µ B /T < ∼ . ≤ √ s NN ≤
200 GeV. Of course, the Taylor expansions will break down, should the elusive criticalpoint in the QCD phase diagram [16, 17] turn out to be present in this range of baryon chemical potentials. Theconvergence of the series thus needs to be monitored carefully.This paper is organized as follows. In the next section we briefly discuss Taylor series for a HRG model in Boltzmannapproximation. This helps to argue for the significance of sixth order Taylor expansions. In Section III, we presentthe basic framework of Taylor series expansions, introduce expansions in the presence of global constraints and discusssome details of our calculations and the ensembles used. In Section IV we discuss the 6 th order Taylor expansion ofQCD thermodynamics in the simplified case of vanishing strangeness and electric charge chemical potentials. Section Vis devoted to the corresponding discussion of strangeness neutral systems n S = 0 with fixed net electric charge ( n Q ) tonet baryon-number ( n B ) ratio, which is of relevance for the description of hot and dense matter formed in heavy ioncollisions where typically n Q /n B (cid:39) .
4. We discuss the relevance of a non-vanishing electric charge chemical potentialby considering electric charge neutral ( n Q /n B = 0) as well as isospin symmetric ( n Q /n B = 1 /
2) systems. At the endof this section we present a parametrization of the equation of state that can easily be used as input for the modelingof the thermal conditions met in heavy ion collisions. In Section VI we present results on lines of constant pressure,energy density and entropy density and compare their dependence on µ B with empirical results for the freeze-outconditions observed in heavy ion collisions. We comment on the radius of convergence of the Taylor series for thepressure and resulting constraints for the location of a possible critical point in Section VII. Finally we present ourconclusions in Section VIII. Details on (A) the statistics and simulation parameters, (B) explicit expressions for theexpansions of electric charge and baryon number chemical potentials, and (C) explicit expressions for the expansionparameters of the lines of constant physics are given in three Appendices A-C. II. TAYLOR EXPANSIONS AND THE LOW AND HIGH TEMPERATURE LIMITS OF STRONGINTERACTION MATTER
The main aim of this work is to supply a EoS of strong-interaction matter using up to sixth order Taylor expansionsfor bulk thermodynamic observables. As we will see later at present results on sixth order expansion coefficients inthe Taylor series will mainly help to constrain truncation errors in the fourth order expansion rather than providingaccurate results on the sixth order contribution to thermodynamic quantities. We will argue that our analysis providesreliable results for the EoS for baryon chemical potentials up to µ B /T (cid:39) T (cid:39)
160 MeV andfor an even larger range in µ B /T at higher temperatures.Before turning to a discussion of lattice QCD results on the EoS, it may be useful to analyze truncation effects inthe hadron resonance gas (HRG) model, which seems to provide a good approximation for thermodynamics in thelow temperature, hadronic regime. For simplicity let us consider the case of vanishing electric charge and strangenesschemical potentials, µ Q = µ S = 0. At temperatures close to the transition temperature T c (cid:39)
154 MeV and for baryonchemical potentials less than a few times the transition temperature, the baryon sector of a HRG is well described inthe Boltzmann approximation. In a HRG model calculation based on non-interacting hadrons the pressure may thenbe written as P ( T, µ B ) = P M ( T ) + P B ( T, ˆ µ B )= P M ( T ) + P B ( T,
0) + P B ( T,
0) (cosh(ˆ µ B ) − , (1)where we introduced the notation ˆ µ B ≡ µ B /T and P M ( T ) ( P B ( T, ˆ µ B )) denote the meson (baryon) contributions tothe pressure. A similar relation holds for the energy density, (cid:15) ( T, µ B ) = (cid:15) M ( T ) + (cid:15) B ( T, ˆ µ B )= (cid:15) M ( T ) + (cid:15) B ( T,
0) + (cid:15) B ( T,
0) (cosh(ˆ µ B ) − , (2)with (cid:15) M/B ≡ T (cid:0) ∂ ( P M/B /T ) /∂T (cid:1) ˆ µ B . The µ B -dependent contribution thus is simple and can easily be representedby a Taylor series. Truncating this expansion at (2 n )-th order we obtain(∆( P/T )) n ≡ ( P B ( T, µ B ) − P B ( T, n T = n (cid:88) k =1 χ B,HRG k ( T )(2 k )! ˆ µ kB (cid:39) P B ( T, T n (cid:88) k =1 k )! ˆ µ kB , (3)where in the last equality we made use of the fact that in HRG models constructed from non-interacting, point-likehadrons, all expansion coefficients are identical when using a Boltzmann approximation for the baryon sector, i.e.all baryon number susceptibilities are identical, χ B,HRG k = P B ( T, χ B,HRG k /χ B,HRG k − = χ B,HRG k /χ B,HRG = 1. Similarly one finds for the net baryon-number density, n B T = P B ( T, T sinh ˆ µ B = ∞ (cid:88) k =1 χ B,HRG k ( T )(2 k − µ k − B (cid:39) P B ( T, T ∞ (cid:88) k =1 k − µ k − B . (4)Higher order corrections are thus more important in the net baryon-number density than in the expansions of thepressure or energy density. For instance, the contribution to µ B n B /T at O (ˆ µ kB ) is a factor 2 k larger than thecorresponding O (ˆ µ kB ) expansion coefficient of the pressure.In Fig. 1 we show results from a Taylor series expansion of the µ B -dependent part of the pressure in a HRG modeltruncated after leading order (LO), next-to-leading order (NLO) and next-to-next-to-leading (NNLO) order. Thesetruncated expansions are compared to the exact result, i.e. (∆ P ) ∞ ( T ) = P B ( T )(cosh(ˆ µ B ) − n -th order truncated Taylor series ((∆ P ) n ( T )) from the exact result ((∆ P ) ∞ ( T )). As can be seenalready the fourth order Taylor series provides a good approximation for the pressure (and energy as well as entropydensity) of a HRG for all µ B ≤ T . At µ B = 2 T the fourth order Taylor series for the µ B -dependent contribution tothe pressure deviates by less than 5% from the exact result. These deviations are, of course, even smaller in the totalpressure which in the temperature range of interest is dominated by the meson contribution. Even at T = 170 MeV,which certainly is already above the range of applicability of HRG models, the baryonic contribution to the pressure(energy density) amounts only to about 20% (30%). A 5% truncation error in the µ B -dependent contribution to thepressure or energy density thus amounts to less than a 2% effect in the total pressure or energy density. Similarestimates hold for the more general case of non-vanishing µ Q and µ S .Of course, the good convergence properties of the Taylor series for the pressure in HRG models also reflect that theradius of convergence of this series is infinite. If there exists a critical point in the QCD phase diagram one cannotexpect to find that the Taylor series is that well behaved. Still the HRG result provides a benchmark also for the QCDcase. If the radius of convergence of the Taylor series for the QCD pressure is finite and, in particular, smaller than µ B (cid:39) T , one should find large deviations in the generalized susceptibilities from the corresponding HRG results.Ratios of susceptibilities have to grow asymptotically like, χ B,QCD k /χ B,QCD k − ∼ k in order to yield a finite radiusof convergence for a Taylor expansion. We will come back to a discussion of this asymptotic behavior after havingdiscussed our sixth order calculation of Taylor expansion coefficients. x= µ B /T( ∆ P) n /P B (T,0) n= ∞ x= µ B /T( ∆ P) n /P B (T,0) n= ∞ x= µ B /T( ∆ P) n /P B (T,0) n= ∞
64 0 1 2 3 4 5 6 7 8 9 10 0 0.5 1 1.5 2 2.5 3 x= µ B /T( ∆ P) n /P B (T,0) n= ∞
642 0 0.1 0.2 0.3 0.4 0.5 0 0.5 1 1.5 2 2.5 3 x= µ B /T1-( ∆ P) n / ( ∆ P) ∞ x= µ B /T1-( ∆ P) n / ( ∆ P) ∞ x= µ B /T1-( ∆ P) n / ( ∆ P) ∞ x= µ B /T1-( ∆ P) n / ( ∆ P) ∞ Figure 1. n -th order Taylor series, (∆ P ) n for (∆ P ) ∞ = P B ( T, x ) −
1) compared to the exact result. The insertionshows the relative error due to truncation of the Taylor series after n -th order. Note that the sixth order result is hardly visiblebehind the exact result. Let us briefly mention also the high temperature limit. At large values of the temperature, the pressure approachesthat of a massless ideal gas of quarks and gluons. In this limit the pressure is just a second order polynomial in ˆ µ f , P ideal T = 8 π
45 + (cid:88) f = u,d,s (cid:20) π
60 + 12 (cid:16) µ f T (cid:17) + 14 π (cid:16) µ f T (cid:17) (cid:21) , (5)In this limit a fourth order Taylor expansion thus provides the exact results for the basic bulk thermodynamicobservables. This also is correct in leading order perturbation theory, i.e. at O ( g ) [18]. III. OUTLINE OF THE CALCULATIONA. Taylor series in baryon number, electric charge and strangeness chemical potentials
Our goal is the calculation of Taylor expansion coefficients for basic bulk thermodynamic observables of strong-interaction matter in terms of chemical potentials µ X for conserved charges ( X = B, Q, S ). We start with theexpansion of the pressure, P , in terms of the dimensionless ratios ˆ µ X ≡ µ X /T , which are the logarithms of fugacities, PT = 1 V T ln Z ( T, V, ˆ µ u , ˆ µ d , ˆ µ s ) = ∞ (cid:88) i,j,k =0 χ BQSijk i ! j ! k ! ˆ µ iB ˆ µ jQ ˆ µ kS , (6)with χ BQS ≡ P ( T, /T . The chemical potentials for conserved charges are related to the quark chemical potentials( µ u , µ d , µ s ), µ u = 13 µ B + 23 µ Q ,µ d = 13 µ B − µ Q ,µ s = 13 µ B − µ Q − µ S . (7)The expansion coefficients χ BQSijk , i.e. the so-called generalized susceptibilities, can be calculated at vanishing chemicalpotential , χ BQSijk ≡ χ BQSijk ( T ) = ∂P ( T, ˆ µ ) /T ∂ ˆ µ iB ∂ ˆ µ jQ ∂ ˆ µ kS (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ˆ µ =0 . (8)From Eq. 6 it is straightforward to obtain the Taylor series for the number densities, n X T = ∂P/T ∂ ˆ µ X , X = B, Q, S . (9)This only requires knowledge of the expansion coefficients entering the series for
P/T . The energy ( (cid:15) ) and entropy ( s )densities, on the other hand, also require derivatives of the generalized susceptibilities with respect to temperature,which are the expansion coefficients of the trace anomaly,∆( T, ˆ µ B , ˆ µ Q , ˆ µ S ) ≡ (cid:15) − PT = T ∂P/T ∂T = ∞ (cid:88) i,j,k =0 Ξ BQSijk i ! j ! k ! ˆ µ iB ˆ µ jQ ˆ µ kS , (10)with i + j + k even and Ξ BQSijk ( T ) = T d χ BQSijk ( T )d T . (11)With this one finds for the Taylor expansions of the energy and entropy densities, (cid:15)T = ∞ (cid:88) i,j,k =0 Ξ BQSijk + 3 χ BQSijk i ! j ! k ! ˆ µ iB ˆ µ jQ ˆ µ kS , (12) sT = (cid:15) + p − µ B n B − µ Q n Q − µ S n S T = ∞ (cid:88) i,j,k =0 Ξ BQSijk + (4 − i − j − k ) χ BQSijk i ! j ! k ! ˆ µ iB ˆ µ jQ ˆ µ kS . (13) B. Constrained series expansions
In our calculations we generated all generalized susceptibilities up to 6 th order, which are needed to set up the generalTaylor series in terms of the three conserved charge chemical potentials as discussed in the previous subsection. Inthe following we will, however, consider only thermodynamic systems, in which the electric charge and strangenesschemical potentials are fixed by additional constraints and become functions of the baryon chemical potential andtemperature. We only consider constraints that can be fulfilled order by order in the Taylor series expansion. Thatis, for the construction of the 6 th order Taylor series of the pressure in terms of ˆ µ B we need to know the expansionof ˆ µ Q ( T, µ B ) and ˆ µ S ( T, µ B ) up to fifth order in ˆ µ B ,ˆ µ Q ( T, µ B ) = q ( T )ˆ µ B + q ( T )ˆ µ B + q ( T )ˆ µ B + . . . , ˆ µ S ( T, µ B ) = s ( T )ˆ µ B + s ( T )ˆ µ B + s ( T )ˆ µ B + . . . . (14)The above parametrization includes the cases of vanishing electric charge and strangeness chemical potentials, µ Q = µ S = 0, which we are going to discuss in the next section as well as the strangeness neutral case with fixed electriccharge to baryon-number ratio, which we will analyze in Section V. We often suppress the argument ( T ) of the generalized susceptibilities. We also suppress superscripts and subscripts of χ BQSijk wheneverone of the subscripts vanishes, e.g. χ BQSi k ≡ χ BSik . Implementing the constraints specified in Eq. 14 in the Taylor series for the pressure and net conserved-chargenumber densities one obtains series in terms of the baryon chemical potential only, P ( T, µ B ) T − P ( T, T = ∞ (cid:88) k =1 P k ( T )ˆ µ kB , (15) n X T = ∞ (cid:88) k =1 N X k − ˆ µ k − B , X = B, Q, S . (16)Using ˆ µ B d P/T dˆ µ B = ˆ µ B n B T + ˆ µ B dˆ µ Q dˆ µ B n Q T + ˆ µ B dˆ µ S dˆ µ B n S T , (17)and the series expansions of ˆ µ Q and ˆ µ S given in Eq. 14 one easily finds the relation between the expansion coefficientsfor the pressure and number densities, P n = 12 n (cid:32) N B n − + n (cid:88) k =1 (2 k − (cid:16) s k − N S n − k +1 + q k − N Q n − k +1 (cid:17)(cid:33) . (18)When imposing constraints on the electric charge and strangeness chemical potentials, these generally becometemperature dependent functions as indicated in Eq. 14. The temperature derivative of P/T at fixed ˆ µ B in theconstraint case and the partial derivative of P/T at fixed (ˆ µ B , ˆ µ Q , ˆ µ S ), which defines the trace anomaly ∆ (Eq. 10),thus are related through T d P/T d T = ∆ + T ˆ µ (cid:48) Q n Q T + T ˆ µ (cid:48) S n S T , (19)where the (total) temperature derivative d / d T is taken at fixed ˆ µ B and ˆ µ (cid:48) X = dˆ µ X / d T . With this we obtain theTaylor series for the trace anomaly,∆( T, ˆ µ B ) = (cid:15) − PT = (cid:18) (cid:15) − PT (cid:19) ˆ µ B =0 + ∞ (cid:88) n =1 ( T P (cid:48) n ( T ) − h n ( T ))ˆ µ nB , (20)with h n = n (cid:88) k =1 (cid:16) s (cid:48) k − N S n − k +1 + q (cid:48) k − N Q n − k +1 (cid:17) . (21)We also introduce t n = n (cid:88) k =1 (cid:16) s k − N S n − k +1 + q k − N Q n − k +1 (cid:17) . (22)With this the Taylor series expansion of the energy and entropy densities for constraint cases, in which ˆ µ Q and ˆ µ S satisfy Eq. 14, becomes (cid:15) ( T, µ B ) T − (cid:15) ( T, T = ∞ (cid:88) n =1 (cid:15) n ( T )ˆ µ nB , (23) s ( T, µ B ) T − s ( T, T = ∞ (cid:88) n =1 σ n ( T )ˆ µ nB . (24)with (cid:15) n ( T ) = 3 P n ( T ) + T P (cid:48) n ( T ) − h n ( T ) and σ k ( t ) = 4 P n ( T ) + T P (cid:48) n ( T ) − N B n − ( T ) − h n ( T ) − t n ( T ). C. Numerical calculation of generalized susceptibilities up to O ( µ ) The generalized susceptibilities χ BQSijk have been calculated on gauge field configurations generated for (2+1)-flavorQCD using the Highly Improved Staggered Quark (HISQ) action [19] and the tree-level improved Symanzik gaugeaction.All calculations are performed using a strange quark mass m s tuned to its physical value. We performed calculationswith two different light to strange quark mass ratios, m l /m s = 1 /
27 and 1 /
20. The former corresponds to a pseudo-scalar Goldstone mass, which in the continuum limit yields a pion mass m π (cid:39)
140 MeV, the latter leads to a pionmass m π (cid:39)
160 MeV. These parameters are fixed using the line of constant physics determined by HotQCD from the f K scale. Using f K = 155 . / √ a ( β ) at a given value of the gaugecoupling β and the corresponding set of quark masses ( m l , m s ), which in turn fixes the temperature on a lattice withtemporal extent N τ , i.e. T = ( N τ a ) − . More details on the scale determination are given in [4].All calculations have been performed on lattices of size N σ N τ with an aspect ratio N σ /N τ = 4. We performcalculations in the temperature interval T ∈ [135 MeV ,
330 MeV] using lattices with temporal extent N τ = 6 , , T ≤
175 MeV all calculations have been performed with the lighter, physical quark mass ratio m l /m s = 1 /
27. Inthe high temperature region quark mass effects are small and we based our calculations on existing data sets for m l /m s = 1 /
20, which have previously been generated by the HotQCD collaboration and used for the calculation ofsecond order susceptibilities [20]. These data sets have been extended for the calculation of higher order susceptibilities.Gauge field configurations are stored after every 10 th molecular dynamics trajectory of unit length.All calculations of 4 th and 6 th order expansion coefficients have been performed on lattices with temporal extent N τ = 6 and 8. In these cases we gathered a large amount of statistics. At low temperatures we have generated up to1 . N τ = 6 and up to 1 . N τ = 8. At high temperature less than a tenthof this statistics turned out to be sufficient. The 2 nd order expansion coefficients have been calculated on latticeswith four different temporal extends, N τ = 6 , , ,
16. At fixed temperature this corresponds to four differentvalues of the lattice cut-off, which we used to extract continuum extrapolated results for the second order expansioncoefficients. We also extrapolated results for the higher order expansion coefficients to the continuum limit. However,having at hand results from only two lattice spacings for these expansion coefficients we consider these extrapolationsas estimates of the results in the continuum limit.On each configuration the traces of all operators needed to construct up to sixth order Taylor expansion coefficientshave been calculated stochastically. For the calculation of 2 nd and 4 th order expansion coefficients we follow thestandard approach of introducing a non-zero chemical potential in the QCD Lagrangian as an exponential prefactorfor time-like gauge field variables [21], i.e. the chemical potential µ f for quark flavor f is introduced through afactor e µ f a (e − µ f a ) on time-like links directed in the forward (backward) direction. This insures that all observablescalculated are free of ultra-violet divergences. For the calculation of all 6 th order expansion coefficients we use the so-called linear- µ approach [22, 23]. This becomes possible as no ultra-violet divergences appear in 6 th order cumulantsand above. In the linear- µ formulation the number of operators that contribute to cumulants is drastically reducedand their structure is simplified. All operators appearing in the exponential formulation, that involve second or higherorder derivatives of the fermion matrix [11], vanish. The remaining operators are identical in both formulations. Onethus only has to calculate traces of observables that are of the form,Tr M − f M (cid:48) f M − f M (cid:48) f ....M − f M (cid:48) f , where M f is the staggered fermion matrix for light ( f = l ) or strange ( f = s ) quarks, respectively, and M (cid:48) f denotes itsderivative with respect to the flavor chemical potential ˆ µ f . The final error on these traces depends on the noise due tothe use of stochastic estimators for the inversion of the fermion matrices M f , as well as on the gauge noise resultingfrom a finite set of gauge configurations that get analyzed. We analyzed the signal to noise ratio for all traces ofoperators that we calculate and identified the operator D = M − f M (cid:48) f , as being particularly sensitive to the stochasticnoise contribution. This operator has been measured using 2000 random noise vectors. For the calculation of traces ofall other operators we used 500 random noise vectors. We checked that this suffices to reduce the stochastic noise wellbelow the gauge noise. The simulation parameters and the statistics accumulated in this calculation are summarizedin the tables of Appendix A.All fits and continuum extrapolations shown in the following are based on spline interpolations with coefficientsthat are allowed to depend quadratically on the inverse temporal lattice size. Our fitting ansatz and the strategyfollowed to arrive at continuum extrapolated results are described in detail in Ref. [3]. For the current analysis wefound it sufficient to use spline interpolations with quartic polynomials and 3 knots whose location is allowed to varyin the fit range. T [MeV] χ free quark gas T c =(154 +/-9) MeV m s /m l =20 (open)27 (filled)PDG-HRGcont. extrap.N τ =161286 T [MeV] cont. extrap.m s /m l =27, N τ =681216PDG-HRGQM-HRG χ Figure 2. The leading order ( O ( µ B )) correction to the pressure calculated at zero baryon chemical potential. The left handfigure shows the leading order correction in a large temperature range. The right hand part of the figure shows an enlarged viewinto the low temperature region. In addition to the continuum extrapolation of the lattice QCD results we also show resultsfrom HRG model calculations based on all hadron resonances listed by the particle data group (PDG-HRG) and obtained inquark model calculations (QM-PDG). IV. EQUATION OF STATE FOR µ Q = µ S = 0 Let us first discuss the Taylor expansion for bulk thermodynamic observables in the case of vanishing electric chargeand strangeness chemical potentials. This greatly simplifies the discussion and yet incorporates all the features of themore general case. Also the discussion of truncation errors presented in this section carries over to the more generalsituation.
A. Pressure and net baryon-number density
For µ Q = µ S = 0 the Taylor expansion coefficients P n and N B n − , introduced in Eqs. 15 and 16, are simply relatedby P n = 12 n N B n − = 1(2 n )! χ B n . (25)The series for the pressure and net baryon-number density simplify to, P ( T, µ B ) − P ( T, T = ∞ (cid:88) n =1 χ B n ( T )(2 n )! (cid:16) µ B T (cid:17) n = 12 χ B ( T )ˆ µ B (cid:18) χ B ( T ) χ B ( T ) ˆ µ B + 1360 χ B ( T ) χ B ( T ) ˆ µ B + ... (cid:19) , (26) n B T = ∞ (cid:88) n =1 χ B n ( T )(2 n − µ n − B = χ B ( T )ˆ µ B (cid:18) χ B ( T ) χ B ( T ) ˆ µ B + 1120 χ B ( T ) χ B ( T ) ˆ µ B + ... (cid:19) . (27)In Eqs. 26 and 27 we have factored out the leading order (LO) µ B -dependent part in the series for the pressure aswell as the net baryon-number density. This helps to develop a feeling for the importance of higher order contributionsand, in particular, the approach to the HRG limit at low temperatures. Note that all ratios χ B n /χ B are unity in aHRG and, in the infinite temperature, ideal quark gas limit, χ B /χ B = 2 / (3 π ) (cid:39) .
068 is the only non-vanishinghigher order expansion coefficient. From Eqs. 26 and 27 it is evident that contributions from higher order expansioncoefficients become more important in the number density than in the pressure. Relative to the LO result, thecontributions of the NLO and NNLO expansion coefficients for n B /T are a factor two and three larger respectivelythan for the corresponding expansion coefficients in the pressure series.We show the leading order coefficient χ B ( T ) in Fig. 2 and the NLO ( χ B ) and NNLO ( χ B ) coefficients divided by χ B ( T ) in Fig. 3. The left hand part of Fig. 2 shows the leading order contribution χ B in the entire temperatureinterval used in the current analysis. For the LO expansion coefficients we also used data from simulations on 48 × m l /m s = 1 /
20 [3] and generated new ensembles for m l /m s = 1 /
27 at nine
HRG free quark gas m s /m l =20 (open)27 (filled) χ B / χ B T [MeV] cont. est.N τ =68 -2-1 0 1 2 3 140 160 180 200 220 240 260 280 m s /m l =20 (open)27 (filled) χ B / χ B T [MeV] cont. est.N τ =68 -2-1 0 1 2 3 140 160 180 200 220 240 260 280 HRG
Figure 3.
Left:
The ratio of fourth and second order cumulants of net-baryon number fluctuations ( χ B /χ B ) versus temperature. Right: same as the left hand side, but for the ratio of sixth and second order cumulants of net-baryon number fluctuations( χ B /χ B ). The boxes indicate the transition region, T c = (154 ±
9) MeV. Grey bands show continuum estimate. temperature values below T = 175 MeV. Furthermore, we used data on 64 ×
16 lattices at a corresponding set oflow temperature values. These data are taken from an ongoing calculation of higher order susceptibilities performedby the HotQCD Collaboration . This allowed us to update the continuum extrapolation for χ B given in [20]. Thenew continuum extrapolation shown in Fig. 2 is consistent with our earlier results, but has significantly smaller errorsin the low temperature region. In the right hand part of this figure we compare the continuum extrapolated latticeQCD data for χ B with HRG model calculations. It is obvious that the continuum extrapolated QCD results overshootresults obtained from a conventional, non-interacting HRG model calculations with resonances taken from the particledata tables (PDG-HRG) and treated as point-like excitations. We therefore compare the QCD results also with aHRG model that includes additional strange baryons,which are not listed in the PDG but are predicted in quarkmodels and lattice QCD calculations. We successfully used such an extended HRG model (QM-HRG) in previouscalculations [5, 6]. As can be seen in Fig. 2 (left), continuum extrapolated results for χ B agree well with QM-HRGcalculations.As can be seen in the left hand part of Fig. 3, the ratio χ B /χ B approaches unity with decreasing temperature, but issmall at high temperatures where the leading order correction is large. The relative contribution of the NLO correctionthus is largest in the hadronic phase, where χ B /χ B (cid:39)
1. For temperatures
T < ∼
155 MeV we find χ B /χ B ≤ .
8. Therelative contribution of the NLO correction to the µ B -dependent part of the pressure (number density) in the crossoverregion and below thus is about 8% (16%) at µ B /T = 1 and rises to about 33% (66%) at µ B /T = 2. At temperatureslarger than 180 MeV the relative contribution of the NLO correction to pressure and number density at µ B /T = 2 isless than 8% and 16%, respectively.The relative contribution of the O (ˆ µ B ) correction, χ B /χ B , is shown in the right hand part of Fig. 3. The ideal gaslimit for this ratio vanishes. Obviously the ratio is already small for all temperatures T >
180 MeV, i.e. χ B /χ B ≤ . µ B = 2 the correction to the leading order result is less than 2.2% for the µ B -dependent part ofthe pressure and less than 7% for the net baryon-number density. At lower temperatures the statistical errors oncurrent results for χ B /χ B are still large. However, a crude estimate for the magnitude of this ratio at all temperatureslarger than 130 MeV suggests, (cid:12)(cid:12) χ B /χ B (cid:12)(cid:12) ≤
3. In the low temperature, hadronic regime and for ˆ µ B = 2 the O (ˆ µ B )corrections to the µ B -dependent part of the pressure can be about 13%. However, in the total pressure, which alsoreceives large contributions from the meson sector, this will result only in an error of less than 3%. In the calculationof the net baryon-number density, on the other hand, the current uncertainty on O (ˆ µ B ) expansion coefficients resultsin errors of about 40% at temperatures below T (cid:39)
155 MeV. In fact, as discussed already in section II, higher ordercorrections are larger in the Taylor expansion of the number density. From Eq. 25 it follows for the ratio of NLOand LO expansion coefficients, N B /N B = 3 P /P . Clearly better statistics is needed in the low temperature rangeto control higher order corrections to n B /T .In Fig. 4 we show results for the µ B -dependent part of the pressure (left) and the net baryon-number density (right)calculated from Taylor series up to and including LO, NLO and NNLO contributions, respectively. This suggests that We thank the HotQCD Collaboration for providing access to the second order quark number susceptibilities. µ B (cid:39) T results for the pressure at low temperature are well described by a Taylor series truncated at NNLO,while at higher temperature NNLO corrections are small even at µ B (cid:39) T . This also is the case for n B /T , althoughthe NNLO correction is large at low temperatures and, at present, does not allow for a detailed quantitative analysisof the baryon-number density in this temperature range.It also is obvious that the Taylor series for the pressure and n B /T in the temperature range up to T (cid:39)
180 MeVare sensitive to the negative contributions of the 6 th order expansion coefficient. The occurrence of a dip in the sixthorder expansion coefficient of the pressure has been expected to show up on the basis of general scaling arguments forhigher order derivatives of the QCD pressure in the vicinity of the chiral phase transition [24]. It may, however, alsoreflect the influence of a singularity on the imaginary chemical potential axis [25] (Roberge-Weiss critical point [26])on Taylor series of bulk thermodynamic observables in QCD. Even with improved statistics it thus is expected thatthe wiggles, that start to show up in the expansion of pressure and net baryon-number density above µ B /T (cid:39) χ B /χ B at T (cid:39)
160 MeV under control in future calculations thus is of importance for the understandingof this non-perturbative regime of the QCD equation of state in the high temperature phase close to the transitionregion. This also indicates that higher order corrections need to be calculated in order to control the equation of statein this temperature regime. µ B /T=2 µ B /T=2.5 µ B /T=1HRG µ Q = µ S =0 [ P ( T , µ B )- P ( T , ) ]/ T T [MeV]
O(µ B6 )O(µ B4 )O(µ B2 ) µ B /T=2 µ B /T=2.5 µ B /T=1HRG µ Q = µ S =0 n B ( T , µ B ) / T T [MeV]
O(µ B6 )O(µ B4 )O(µ B2 ) Figure 4. The µ B -dependent contribution to the pressure (left) and the baryon-number density (right) in the case of vanishingelectric charge and strangeness chemicals potential for several values of the baryon chemical potential in units of temperature.The different bands show results including Taylor series results upto the order indicated. B. Net strangeness and net electric charge densities
For vanishing strangeness and electric charge chemical potentials the corresponding net strangeness ( n S ) and netelectric charge ( n Q ) densities are nonetheless non-zero because the carriers of these quantum numbers also carrybaryon number. The ratios of number densities are given by n X n B = χ BX + χ BX ˆ µ B + χ BX ˆ µ B χ B + χ B ˆ µ B + χ B ˆ µ B , X = Q, S . (28)In a hadron resonance gas the ratios n S /n B and n Q /n B are independent of the baryon chemical potential and,irrespective of the value of ˆ µ B , these ratios approach − T → ∞ limit. One thus mayexpect that these ratios only show a mild dependence on ˆ µ B , which indeed is apparent from the results of the NNLOexpansions shown in Fig. 5.For µ Q = µ S = 0 non-vanishing electric charge and strangeness densities only arise due to a non-zero baryon-chemical potential. In the low temperature HRG phase n Q and n S thus only receive contributions from chargedbaryons or strange baryons, respectively. The ratios n Q /n B and n S /n B thus are sensitive to the particle content ina hadron resonance gas and a comparison with PDG-HRG and QM-HRG is particularly sensitive to the differencesin the baryon content in these two models. It is apparent from Fig. 5 that at low temperatures the QM-HRG modelprovides a better description of the lattice QCD results than the PDG-HRG model.1 µ Q = µ S =0 - n S ( T , µ B ) / n B ( T , µ B ) T [MeV] µ B /T=2.52.01.0O(µ B0 ) PDG-HRGQM-HRG µ Q = µ S =0 n Q ( T , µ B ) / n B ( T , µ B ) T [MeV] µ B /T=2.52.01.0O(µ B0 )PDG-HRGQM-HRG Figure 5. The ratio of net strangeness and net baryon-number densities (left) and the ratio of net electric charge and netbaryon-number densities (right). At low temperatures results from hadron resonance gas calculations at µ B = 0 are shown (seetext). T [MeV] µ Q = µ S =03 P /T T [MeV] µ Q = µ S =03 P /T ε /T T [MeV] µ Q = µ S =03 P /T ε /T σ /T T [MeV] µ Q = µ S =03 P /T ε /T σ /T -0.02-0.01 0 0.01 0.02 0.03 0.04 140 160 180 200 220 240 260 280 T [MeV] -0.02-0.01 0 0.01 0.02 0.03 0.04 140 160 180 200 220 240 260 280
T [MeV] /T -0.02-0.01 0 0.01 0.02 0.03 0.04 140 160 180 200 220 240 260 280 T [MeV] /T ε /T -0.02-0.01 0 0.01 0.02 0.03 0.04 140 160 180 200 220 240 260 280 T [MeV] /T ε /T σ /T -0.02-0.01 0 0.01 0.02 0.03 0.04 140 160 180 200 220 240 260 280 T [MeV] /T ε /T σ /T -0.02 0 0.02 0.04 0.06 0.08120 140 160 180 200 220 240 260 280ratio of NLO and LOexpansion coe ffi cients -0.02 0 0.02 0.04 0.06 0.08120 140 160 180 200 220 240 260 280ratio of NLO and LOexpansion coe ffi cients -0.02 0 0.02 0.04 0.06 0.08120 140 160 180 200 220 240 260 280ratio of NLO and LOexpansion coe ffi cients -0.02 0 0.02 0.04 0.06 0.08120 140 160 180 200 220 240 260 280ratio of NLO and LOexpansion coe ffi cients -0.02 0 0.02 0.04 0.06 0.08120 140 160 180 200 220 240 260 280ratio of NLO and LOexpansion coe ffi cients -0.02 0 0.02 0.04 0.06 0.08120 140 160 180 200 220 240 260 280ratio of NLO and LOexpansion coe ffi cients -0.02 0 0.02 0.04 0.06 0.08120 140 160 180 200 220 240 260 280ratio of NLO and LOexpansion coe ffi cients Figure 6. Leading order (left) and next-to-leading order (right) expansion coefficients for the µ B -dependent part of pressure,the energy and entropy densities in the case µ Q = µ S = 0. The inset in the right hand figure shows the ratios of NLO andLO expansion coefficients P /P , (cid:15) /(cid:15) and σ /σ . Note that the expansion coefficients for the net baryon-number density aredirectly proportional to those of the pressure series, i.e. N B = 2 P and N B = 4 P . C. The energy and entropy densities
In order to calculate the energy and entropy densities, defined in Eqs. 23 and 24, we need to extract the temperaturederivative of the expansion coefficients of the pressure. We use as a starting point the representation of the pressuregiven in Eq. 26 and calculate the temperature derivatives of χ Bn from the splines used to fit this observable. With thiswe construct the expansion coefficients (cid:15) Bn ( T ) and σ Bn defined in Eqs. 12 and 13,∆ (cid:0) (cid:15)/T (cid:1) = (cid:15) ( T, µ B ) − (cid:15) ( T, T = (cid:88) k =1 (cid:15) k ˆ µ kB = (cid:88) k =1 ( T P (cid:48) k + 3 P k ) ˆ µ kB , (29)∆ (cid:0) s/T (cid:1) = s ( T, µ B ) − s ( T, T = (cid:88) k =1 σ k ˆ µ kB = (cid:88) k =1 ( (cid:15) k − (2 k − P k )ˆ µ kB . (30)We show the LO and NLO expansion coefficients for energy and entropy densities together with the expansioncoefficient for the pressure in Fig. 6. Because of Eq. 25 the expansion coefficients of the net baryon-number densityare simply proportional to those of the pressure.Clearly the temperature dependence of the expansion coefficients of the energy and entropy densities shows more2 T [MeV]P(T, µ B )/T µ Q = µ S =0 µ B /T=2.5 0 0.5 1 1.5 2 2.5 3 3.5 4 140 160 180 200 220 240 260 280 T [MeV]P(T, µ B )/T µ Q = µ S =0 µ B /T=2.52.0 0 0.5 1 1.5 2 2.5 3 3.5 4 140 160 180 200 220 240 260 280 T [MeV]P(T, µ B )/T µ Q = µ S =0 µ B /T=2.52.01.0 0 0.5 1 1.5 2 2.5 3 3.5 4 140 160 180 200 220 240 260 280 T [MeV]P(T, µ B )/T µ Q = µ S =0 µ B /T=2.52.01.00.0 0 2 4 6 8 10 12 14 16 140 160 180 200 220 240 260 280 T [MeV] ε (T, µ B )/T µ Q = µ S =0 µ B /T=2.5 0 2 4 6 8 10 12 14 16 140 160 180 200 220 240 260 280 T [MeV] ε (T, µ B )/T µ Q = µ S =0 µ B /T=2.52.0 0 2 4 6 8 10 12 14 16 140 160 180 200 220 240 260 280 T [MeV] ε (T, µ B )/T µ Q = µ S =0 µ B /T=2.52.01.0 0 2 4 6 8 10 12 14 16 140 160 180 200 220 240 260 280 T [MeV] ε (T, µ B )/T µ Q = µ S =0 µ B /T=2.52.01.00.0 Figure 7. (Left) The total pressure in (2+1)-flavor QCD in O (ˆ µ B ) for several values of µ B /T . (Right) The total energy densityin (2+1)-flavor QCD in O (ˆ µ B ) for several values of µ B /T . The results for ˆ µ B = 0 are taken from Ref. [3]. structure than in the case of the pressure. Qualitatively this can be understood in terms of the pseudo-criticalbehavior of bulk thermodynamic observables. Once thermodynamic quantities are dominated by contributions fromthe singular part of the free energy, which is expected to happen in the transition region, they become functions of( T − T c ) + κ ˆ µ B . The temperature derivative of the expansion coefficient P , which gives (cid:15) , thus will show propertiessimilar to those of P . The LO correction (cid:15) B /T has a mild peak, which results from the strongly peaked T -derivativeof χ B which is qualitatively similar to χ B , and the NLO correction is negative in a small temperature interval above T c , which arises from the negative T -derivative of χ B at high temperature, which resembles the negative part of χ B at high temperature.Although the temperature dependence of (cid:15) n and σ n differs from that of the pressure coefficient, P n , the conclusionsdrawn for the relative strength of the expansion coefficients are identical in all cases. As can be seen from the insetin Fig. 6 (right) the relative contribution of the NLO expansion coefficients never exceeds 10%. In particular, attemperatures larger than 180 MeV the magnitude of the NLO expansion coefficients never exceeds 2% of the LOexpansion coefficients. Again this leads to the conclusion that at µ B /T = 2 and temperatures above 180 MeV theNLO correction contributes less than 8% of the leading correction to µ B -dependent part of the energy and entropydensities. For T < ∼
155 MeV, however, the NLO contribution can rise to about 30%. A similar conclusion holds forthe O (ˆ µ B ) corrections, although it requires higher statistics to better quantify the magnitude of this contribution. InFig. 7 we show results for the total pressure and total energy density. For P/T and (cid:15)/T at µ B = 0 we used theresults obtained by the HotQCD Collaboration [3] and added to it the results from the O (ˆ µ B ) expansions presentedabove. This figure also makes it clear that despite of the large error of higher order expansion coefficients, which wehave discused above, the error on the total pressure and energy density still is dominated by errors on their values at µ B = 0. V. EQUATION OF STATE IN STRANGENESS NEUTRAL SYSTEMSA. Taylor expansion of pressure, baryon-number, energy and entropy densities
We now want to discuss the equation of state for strangeness neutral systems with a fixed ratio of electric chargeto baryon-number density, i.e. we impose the constraints [27] n S = 0 , n Q n B = r . (31)These constraints can be realized through suitable choices of the electric charge and strangeness chemical potentials.This thus is a particular case of the constraint expansion discussed in Subsection III B. The expansion coefficients q n , s n , n = 1 , , r = 0 . r , including the case of isospin symmetric systems ( r = 1 /
2) andelectric charge neutral matter ( r = 0).3 P T[MeV]
HRGcont. extrap.N τ =161286 m s /m l =20 (open)27 (filled)free quark gasn S =0, n Q /n B =0.4 P T[MeV]
HRGcont. est.N τ = 86 m s /m l =20 (open)27 (filled)free quark gasn S =0, n Q /n B =0.4 -0.00015-0.0001-5x10 -5 -5 P T[MeV]
HRGcont. est.N τ = 86 -0.00015-0.0001-5x10 -5 -5 m s /m l =20 (open)27 (filled)n S =0, n Q /n B =0.4 S =0, n Q /n B =0.4 N B n - / n P n T [MeV] N B1 /2P cont. extrap.N τ =1286 Figure 8. Expansion coefficients of the pressure (top, and bottom left) and the ratio of net baryon-number density and pressureexpansion coefficients (bottom, right) in strangeness neutral systems with r = 0 .
4. Broad bands show continuum extrapolationsas discussed in Section III. The darker lines in the center of the error bands of these extrapolations show the interpolating fitsdiscussed in Subsection V B. At low temperature lines for HRG model calculations based on hadron resonances listed by theParticle Data Group is shown.
Using the constraints specified in Eq. 31 and the definition of the pressure in terms of generalized susceptibilities, χ BQSijk , the expansion coefficients P n can easily be determined. Here it advantageous to use the relation betweenthe Taylor expansion coefficients of the pressure, P n , and number densities, N X n − , given in Eq. 18, which simplifiesconsiderably for strangeness neutral systems. It now involves only the net baryon-number density coefficients, P = 12 (cid:2) N B + rq N B (cid:3) , (32) P = 14 (cid:2) N B + r (cid:0) q N B + 3 q N B (cid:1)(cid:3) , (33) P = 16 (cid:2) N B + r (cid:0) q N B + 3 q N B + 5 q N B (cid:1)(cid:3) . (34)Explicit expressions for all N Bn − and q n − , for n = 2 , ,
6, are given in Appendix B. The resulting expansioncoefficients for the pressure are shown in Fig. 8. Also shown in the bottom-right panel of this figure is the ratio of theexpansion coefficients for the net baryon-number density, N Bn − and the appropriately rescaled expansion coefficientsof the pressure, nP n . In electric charge neutral systems, r = 0 as well as in the isospin symmetric limit r = 1 /
2, forwhich the expansion coefficients q i = 0 vanish for all i , this ratio is unity. In both cases the simple relation given inEq. 25 holds. Also for other values of r the contribution from terms proportional to r are small. In Fig. 8 (bottom,right) we show the ratio N B n − /nP n for the case r = 0 . n = 2 and 4, respectively. At O (ˆ µ B ) differences between N B and 2 P never exceed 2% and at O (ˆ µ B ) the difference between N B and 4 P varies between 3% at low temperatureand -6% at high temperature. In the infinite temperature ideal gas limit the ratios become N B / P = 1 .
018 and N B / P = 0 . µ S =0, µ Q =0n S =0, r=0.00.20.4 T[MeV]P /P m s /m l =20 (open)27 (filled) Figure 9. Ratio of O (ˆ µ B ) expansion coefficients of the pressure in systems with electric charge to net baryon-number ratio r = n Q /n B relative to that of strangeness neutral, isospin symmetric systems ( r = 1 / N τ = 8. number-ratio is weak. The O (ˆ µ B ) expansion coefficient of the pressure in strangeness neutral systems differs by atmost 10% in electric charge neutral ( r = 0) and isospin symmetric systems ( r = 1 / P evaluated for different values of r is shown in Fig. 9. For chemical potentials ˆ µ ≤ O (ˆ µ B ) expansion coefficients differ byalmost 50% in the high temperature limit. For T <
150 MeV this difference is only about 10% reflecting that thedifferent treatment of the strangeness sector becomes less important for the thermodynamics at low temperature.This is also shown in Fig. 9. P / P T[MeV]
HRGcont. est.N τ = 86 m s /m l =20 (open)27 (filled)free quark gasn S =0, n Q /n B =0.4 -2-1 0 1 2 140 160 180 200 220 240 260 280 P / P T[MeV]
HRGcont. est.N τ = 86 -2-1 0 1 2 140 160 180 200 220 240 260 280 m s /m l =20 (open)27 (filled)n S =0, n Q /n B =0.4 Figure 10. Ratio of expansion coefficients of the pressure in strangeness neutral systems with r = 0 .
4. The darker lines in thecenter of the error bands of these extrapolations show results obtained with the parametrization discussed in Subsection V B.
Compared to the leading O (ˆ µ B ) contributions to bulk thermodynamic observables the O (ˆ µ B ) and O (ˆ µ B ) correctionsare smaller in the strangeness neutral case than in the case µ Q = µ S = 0, which we have discussed in the previoussection. This is evident from Fig. 10, where we show the ratios 12 P /P and 360 P /P . These combinations areunity in a HRG with µ S = µ Q = 0 but smaller than unity in the strangeness neutral case. Higher order correctionsin Taylor series for strangeness neutral systems thus are of less importance than in the case µ S = 0. This alsomeans that the errors, which are large on e.g. sixth order expansion coefficients, are of less importance for the overallerror budget of Taylor expansions in strangeness neutral systems. This is indeed reflected in the µ B -dependence of5 µ B /T=2 µ B /T=2.5 µ B /T=1HRGn S =0, n Q /n B =0.4 [ P ( T , µ B )- P ( T , ) ]/ T T [MeV] O( µ B6 )O( µ B4 )O( µ B2 ) µ B /T=2 µ B /T=2.5 µ B /T=1HRG n S =0, n Q /n B =0.4 n B ( T , µ B ) / T T [MeV] O( µ B6 )O( µ B4 )O( µ B2 ) µ B /TP(T, µ B )-P(T,0)P HRG (T, µ B )-P HRG (T,0)T=145 MeV155 MeV165 MeV µ B /Tn B (T, µ B )n BHRG (T, µ B )T=145 MeV155 MeV165 MeV Figure 11. The µ B dependent contribution to the pressure (top, left) and the baryon-number density (top, right) for severalvalues of the baryon chemical potential in units of temperature. The lower two panels show these quantities normalized to thecorresponding HRG model values, obtained from a calculation with all baryon resonances, up to mass m H = 2 . µ B /T for three values of the temperature. ( P ( T, µ B ) − P ( T, /T and n B ( T, µ B ) /T shown in the upper panels of Fig. 11 for the case r = 0 .
4. As can be seenin these two figures, at low temperatures the µ B -dependent part of the pressure as well as the net baryon-numberdensity agree quite well with HRG model calculations that describe the thermodynamics of a gas of non-interacting,point-like hadron resonances. This agreement, however, gets worse at larger values of µ B . Not unexpectedly, athigher temperatures deviations from HRG model calculations become large already at small values of µ B . This isapparent from the lower two panels of Fig. 11, where we show the ratio of the µ B -dependent part of the pressureand the corresponding HRG model result (left) and the net baryon-number density divided by the correspondingHRG model result (right). In the HRG model calculation ( P ( T, µ B ) − P ( T, /T as well as n B ( T, µ B ) /T onlydepend on the baryon sector of the hadron spectrum. The results shown in Fig. 11 thus strongly suggest that HRGmodel calculations using resonance spectra in model calculations for non-interacting, point-like hadron gases may beappropriate (within ∼
10% accuracy) to describe the physics in the crossover region of strongly interacting matterat vanishing or small values of the baryon chemical potential, but fail to do so at large µ B /T and/or T > ∼
160 MeV.At T = 165 MeV QCD and HRG model results for the net baryon-number density differ by 40% at µ B /T = 2. Thishas consequences for the determination of freeze-out conditions in heavy ion collisions. We will come back to thisdiscussion in Section VI.The µ B -dependent contributions to the energy and entropy densities have been defined in Eqs. 23 and 24. In It has been pointed out that the point-like particle approximation is appropriate in the meson sector but not in the baryon sector athigh density. Introducing a non-zero size of hadron resonances [28, 29] may, for some observables, improve the comparison with QCDthermodynamics [30, 31]. However, it seems that the introduction of several additional parameters will be needed to achieve overallgood agreement with the many observables calculated now in QCD in the temperature range of interest, i.e. in the crossover regionfrom a hadron gas to strongly interacting quark-gluon matter. T [MeV]n S =0, n Q /n B =0.4 /T T [MeV]n S =0, n Q /n B =0.4 /T ε /T T [MeV]n S =0, n Q /n B =0.4 /T ε /T σ /T T [MeV]n S =0, n Q /n B =0.4 /T ε /T σ /T T [MeV]n S =0, n Q /n B =0.4 /T ε /T σ /T T [MeV]n S =0, n Q /n B =0.4 /T ε /T σ /T -0.015-0.01-0.005 0 0.005 0.01 0.015 0.02 0.025 140 160 180 200 220 240 260 280 T [MeV] /T -0.015-0.01-0.005 0 0.005 0.01 0.015 0.02 0.025 140 160 180 200 220 240 260 280 T [MeV] /T ε /T -0.015-0.01-0.005 0 0.005 0.01 0.015 0.02 0.025 140 160 180 200 220 240 260 280 T [MeV] /T ε /T σ /T -0.015-0.01-0.005 0 0.005 0.01 0.015 0.02 0.025 140 160 180 200 220 240 260 280 T [MeV] /T ε /T σ /T -0.015-0.01-0.005 0 0.005 0.01 0.015 0.02 0.025 140 160 180 200 220 240 260 280 T [MeV] /T ε /T σ /T -0.015-0.01-0.005 0 0.005 0.01 0.015 0.02 0.025 140 160 180 200 220 240 260 280 T [MeV] /T ε /T σ /T -0.04 0 0.04 0.08 0.12 0.16 140 160 180 200 220 240 260 280ratio of NLO and LOexpansion coe ffi cients Figure 12. Leading order (left) and next-to-leading order (right) expansion coefficients for the µ B -dependent part of pressure,the energy and entropy densities in the strangeness neutral case with fixed electric charge to net baryon-number density, n Q /n B = 0 .
4. The darker lines in the center of the error bands of these extrapolations show the interpolating fits discussed inSubsection V B. The insert in the right hand figure shows the ratios of NLO and LO expansion coefficients N B /N B , P /P , (cid:15) /(cid:15) and σ /σ . The influence of a non-vanishing electric charge chemical potential, which formally gives rise to deviationsfrom the result in the isospin symmetric limit ( N B = 2 P , N B = 4 P ), are negligible at O (ˆ µ B ) and O (ˆ µ B ). For that reasonwe do not show results for N B and N B . However, we show in the insertion in the left hand figure the ratio N B /N B (blackline) which clearly shows that NLO corrections are a factor two larger in the Taylor series for the number density then in thepressure series. strangeness neutral systems the expansion coefficients simplify considerably, (cid:15) n ( T ) = 3 P n ( T ) + T P (cid:48) n ( T ) − r n (cid:88) k =1 T q (cid:48) k − N B n − k +1 (35) σ n ( T ) = 4 P n ( T ) + T P (cid:48) n ( T ) − N B n − − r n (cid:88) k =1 ( q k − + T q (cid:48) k − ) N B n − k +1 (36)Results for the O ( µ B ) and O ( µ B ) expansion coefficients are shown in Fig. 12 together with the corresponding expansioncoefficients for the pressure and net baryon-number density. Results for the total energy density as well as the totalpressure for µ B /T = 0 and 2 are shown in Fig. 13. As discussed in the previous section also here it is evident thatcurrent errors on the total pressure and energy density are dominated by errors on these observables at µ B = 0.In Fig. 13 we also show results for the total pressure obtained within the stout discretization scheme. The result forˆ µ B = 0 is taken from [2]. The ˆ µ B -dependent contribution is based on calculations with an imaginary chemical potential[14]. These results have been analytically continued to real values of ˆ µ B using a 6 th order polynomial in ˆ µ B . As canbe seen the total pressure agrees quite well with the results obtained with a sixth order Taylor expansion, althoughthe results obtained the analytic continuation within the stout discretization scheme tend to stay systematically belowthe central values obtained from the analysis of Taylor series expansions in the HISQ discretization scheme. B. Parametrization of the equation of state At µ B = 0 the HotQCD Collaboration presented a parametrization of the pressure, obtained as interpolating curvesfor the continuum extrapolated fit, that also provided an adequate description of all the other basic thermodynamicquantities, i.e. the energy and entropy densities as well as the specific heat and the velocity of sound [3]. Here wewant to extend this parametrization to the case ˆ µ B >
0. Similar to what has been done at µ B = 0 it turns outthat a ratio of fourth order polynomials in the inverse temperature is flexible enough to describe the temperaturedependence of all required Taylor expansion coefficients in the temperature range T ∈ [130 MeV ,
280 MeV]. We usesuch an ansatz for the three expansion coefficients of the net baryon-number density ( N B , N B , N B ) and the threeelectric charge chemical potentials ( q , q , q ). This suffices to calculate all thermodynamic observables in strangenessneutral systems.We use a ratio of fourth order polynomials in 1 /T as an ansatz for the expansion coefficients of the net baryon-7 S =0 , n Q /n B =0,4 stout, imag. µ B HISQ, real µ B T [MeV] ε /T : µ B /T= 2.003P/T : µ B /T= 2.0 0 Figure 13. The total energy density (upper two curves) of (2+1)-flavor QCD for µ B /T = 0 and 2, respectively. The lower twocurves show corresponding results for three times the pressure. The dark lines show the results obtained with the stout actionfrom analytic continuation with sixth order polynomials in ˆ µ B [14]. number density, N Bk ( T ) = N Bk, n + N Bk, n ¯ t + N Bk, n ¯ t + N Bk, n ¯ t + N Bk, n ¯ t N Bk, d ¯ t + N Bk, d ¯ t + N Bk, d ¯ t + N Bk, d ¯ t , k = 1 , , . (37)Here ¯ t = T c /T and the QCD transition temperature T c = 154 MeV is used as a convenient normalization. Similarlywe define the parametrization of the expansion coefficients for the electric charge chemical potential, q k ( T ) = q k, n + q k, n ¯ t + q k, n ¯ t + q k, n ¯ t + q k, n ¯ t q k, d ¯ t + q k, d ¯ t + q k, d ¯ t + q k, d ¯ t , k = 1 , , . (38)The parameters for these interpolating curves are summarized in Table I.The expansion coefficients of the pressure are then obtained by using Eqs. 32-34. The resulting interpolatingcurves for P k are shown as darker curves in Fig. 8. All other interpolating curves shown as darker curves in otherfigures have been obtained by using the above interpolations. In particular, interpolating curves for the energy andentropy densities are obtained by using Eqs. 35 and 36 and calculating analytically temperature derivatives of theparametrizations of P n and q n given in Eqs. 37 and 38. The resulting interpolating curves for the second and fourthorder Taylor expansion coefficients are shown in Fig. 12.We also used a ratio of fourth order polynomials to interpolate results for the pressure at µ B = 0. We write thepressure as P ( T, µ B = 0) T = p n + p n ¯ t + p n ¯ t + p n ¯ t + p n ¯ t p d ¯ t + p d ¯ t + p d ¯ t + p d ¯ t . (39)The coefficients p in and p id are also given in Table I. VI. LINES OF CONSTANT PHYSICS TO O ( µ B ) We want to use here the Taylor series for bulk thermodynamic observables, i.e. the pressure, energy and entropydensities, to discuss contour lines in the T - µ B plane on which these observables stay constant. It has been argued quitesuccessfully that the thermal conditions at the time of chemical freeze-out in heavy ion collisions can be characterizedby lines in the T - µ B plane on which certain thermodynamic observables or ratios thereof stay constant [32, 33],although the freeze-out mechanism in the rapidly expanding fireball created in a heavy ion collision is of dynamicalorigin and will in detail be more complicated (see for instance [34]). While lines of constant physics (LCPs) involvingtotal baryon-number densities, as used in [32, 33], are not appropriate for calculations within the framework ofquantum field theories, other criteria like lines of constant entropy density in units of T [35] or constant pressure[36–38] have been suggested to characterize freeze-out parameters ( T f , µ fB ) corresponding to heavy ion collisions at8 k N Bk, n N Bk, n N Bk, n N Bk, n N Bk, n N Bk, d N Bk, d N Bk, d N Bk, d k q k, n q k, n q k, n q k, n q k, n q k, d q k, d q k, d q k, d p n p n p n p n p n p d p d p d p d n Q /n B = 0 .
4. These interpolations have been determined for the temperature interval T ∈ [130 MeV ,
280 MeV]. Alsogiven are parameters needed for the interpolation of the expansion coefficients for the electric charge chemical potential (Eq. 38)and the coefficients for the parametrization of the pressure at µ B = 0 given in Eq. 39. different values of the beam energy ( √ s NN ). Generally such criteria have been established by comparing experimentaldata with model calculations based on some version of a HRG model. We will determine here LCPs from the latticeQCD calculations of pressure, energy and entropy densities and confront them with freeze-out parameters that havebeen obtained by comparing particle yields, measured at different values of √ s NN , to HRG model calculations.We consider an observable f ( T, µ B ), i.e. the pressure, energy density or entropy density which are even functionsof µ B . We parametrize a ’line of constant f ’ by, T f ( µ B ) = T (cid:32) − κ f (cid:18) µ B T (cid:19) − κ f (cid:18) µ B T (cid:19) (cid:33) . (40)In order to determine the expansion coefficients κ f and κ f we need to expand the function f ( T, µ B ) up to 4 th orderin µ B and up to second order in T around some point ( T , f ( T, µ B ) = f ( T ,
0) + ∂f ( T, µ B ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) ( T , ( T − T ) + 12 ∂ f ( T, µ B ) ∂µ B (cid:12)(cid:12)(cid:12)(cid:12) ( T , µ B (41)+ 12 ∂ f ( T, µ B ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) ( T , ( T − T ) + 12 ∂∂T ∂ f ( T, µ B ) ∂µ B (cid:12)(cid:12)(cid:12)(cid:12) ( T , ( T − T ) µ B + 14! ∂ f ( T, µ B ) ∂µ B (cid:12)(cid:12)(cid:12)(cid:12) ( T , µ B . Note that we expand here in terms of µ B rather than in ˆ µ B ≡ µ B /T . Replacing the temperature T in Eq. 41 by theansatz for a line of constant f , Eq. 40, and keeping terms up to O ( µ B ) gives f ( T ( µ B ) , µ B ) = f ( T ,
0) + (cid:32) − κ f ∂f ( T, µ B ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) ( T , T + 12 ∂ f ( T, µ B ) ∂µ B (cid:12)(cid:12)(cid:12)(cid:12) ( T , (cid:33) µ B + (cid:32) − κ f ∂f ( T, µ B ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) ( T , T + 12 ∂ f ( T, µ B ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) ( T , ( κ f ) T − ∂∂T ∂ f ( T, µ B ) ∂µ B (cid:12)(cid:12)(cid:12)(cid:12) ( T , κ f T + 14! ∂ f ( T, µ B ) ∂µ B (cid:12)(cid:12)(cid:12)(cid:12) ( T , (cid:33) µ B . We then can determine κ f and κ f by demanding that the expansion coefficients at O ( µ B ) and O ( µ B ) vanish, i.e. κ f = T ∂ f ( T,µ B ) ∂µ B (cid:12)(cid:12)(cid:12) ( T , ∂f ( T,µ B ) ∂T (cid:12)(cid:12)(cid:12) ( T , , (42) κ f = T ∂ f ( T,µ B ) ∂T (cid:12)(cid:12)(cid:12) ( T , ( κ f ) − T ∂∂T ∂ f ( T,µ B ) ∂µ B (cid:12)(cid:12)(cid:12) ( T , κ f + T ∂ f ( T,µ B ) ∂µ B (cid:12)(cid:12)(cid:12) ( T , T ∂f ( T,µ B ) ∂T (cid:12)(cid:12)(cid:12) ( T , . (43)9As we will deal with observables that are given as a Taylor series in ˆ µ B at fixed T , i.e. f ( T, µ B ) = (cid:80) ∞ k =0 f k ( T )ˆ µ kB , thederivatives with respect to µ B appearing in Eqs. 42 and 43 can be replaced by suitable Taylor expansion coefficientsof f ( T, µ B ), κ f = f ( T ) T ∂f ( T ) ∂T (cid:12)(cid:12)(cid:12) ( T , (44) κ f = T ∂ f ( T ) ∂T (cid:12)(cid:12)(cid:12) ( T , ( κ f ) − (cid:16) T ∂f ( T ) ∂T (cid:12)(cid:12)(cid:12) ( T , − f ( T ) (cid:19) κ f + f ( T ) T ∂f ( T ) ∂T (cid:12)(cid:12)(cid:12) ( T , (45)We will in the following work out detailed expressions for the quadratic correction coefficient, κ f , for lines of constantpressure ( f ≡ P ), energy density ( f ≡ (cid:15) ) and entropy density ( f ≡ s ) in strangeness neutral systems with electriccharge to net baryon-number ratio r = 0 .
4. Details for the quartic coefficient, κ f , are given in Appendix C.pressure f ≡ P :The function f ( T, µ B ) is given by P = T (cid:80) n P n ( µ B /T ) n , with P = P ( T, /T denoting the pressure inunits of T at vanishing baryon chemical potential and P n ( T ), n >
0, denoting the expansion coefficients of P ( T, µ B ) /T as introduced in Eq. 15. In the denominator of Eq. 44 we use the thermodynamic relation betweenpressure and entropy density s = ( ∂P/∂T ) µ B . The numerator is given by f ( T ) = T P ( T ). This gives κ P = P s/T , (46)where s/T is evaluated at ˆ µ B = 0.energy density f ≡ (cid:15) :The function f ( T, µ B ) is given by (cid:15) = T (cid:80) n (cid:15) n ( µ B /T ) n , with (cid:15) = (cid:15) ( T, /T denoting the energy density inunits of T at vanishing baryon chemical potential . In the denominator of Eq. 44 we use the thermodynamicrelation between energy density and specific heat C V = ( ∂(cid:15)/∂T ) µ B In the numerator we have f ( T ) = T (cid:15) ( T ).This gives κ (cid:15) = (cid:15) C V /T , (47)where C V /T is evaluated at ˆ µ B = 0.entropy density f ≡ s :The function f ( T, µ B ) is given by s = ( (cid:15) + P − µ B n B − µ Q n Q ) /T = ( (cid:15) + P − µ B n B (1 + rµ Q /µ B )) /T . As n B is of O ( µ B ) we need for the ratio of electric charge and strangeness chemical potentials only the leading orderrelation µ Q /µ B = q defined in Eq. 14. In the denominator we use, ∂s∂T = ∂ ( (cid:15) + P ) /T∂T = − sT + 1 T ∂ ( (cid:15) + P ) ∂T = C V T . (48)In the numerator we have f ( T ) = T ( (cid:15) + P − N B (1 + rq )). With this we get, κ s = T (cid:15) + P − N B (1 + rq ) C V = (cid:15) − P C V /T . (49)where we have used Eq. 32 to replace N B in favor of P .We note that κ (cid:15) > κ s , i.e. with increasing µ B the entropy density decreases on lines of constant energy density.The second order coefficients for the lines of constant physics thus can directly be calculated using the continuumextrapolated results for the pressure and energy density obtained at vanishing chemical potential in [3] and the leadingorder expansion coefficient of the pressure shown in Fig. 10. Similarly we obtain the quartic coefficients from thefourth order expansion of the pressure using the relations given in Appendix C. We show results for κ f and κ f inFig. 14.0 κ κ ε κ HRG0.0000.0020.0040.0060.0080.0100.0120.014130 140 150 160 170 180 190 200T [MeV] κ κ κ ε κ HRG0.00000.00010.0002130 140 150 160 170 180 190 200T [MeV] κ Figure 14.
Left:
Second order curvature coefficients of lines of constant pressure, energy density and entropy density versustemperature in (2+1)-flavor QCD (bands) and in a HRG model (lines).
Right: same as on the left, but for fourth ordercoefficients. The darker lines in the center of the error bands show the interpolating fits discussed in subsection V B. For κ (cid:15) and κ s only these interpolating curves are shown. STAR ALICE Becattini et al. T f ( µ B ) [ M e V ] constant P ε scrossover lines µ B [MeV] µ B /Tn B (T, µ B ) [fm -3 ] ε =556(57) MeV/fm (top) (middle) (bottom) constant P ε s Figure 15.
Left:
Lines of constant pressure, energy density and entropy density versus temperature in (2+1)-flavor QCD forthree different initial sets of values fixed at µ B = 0 and T = 145 MeV, 155 MeV and 165 MeV, respectively (see Table II). Datapoints show freeze-out temperatures determined by the STAR Collaboration in the BES at RHIC (squares) [39] and the ALICECollaboration at the LHC (triangle) [40]. The circles denote hadronization temperatures obtained by comparing experimentaldata on particle yields with a hadronization model calculation [41]. Also shown are two lines representing the current spreadin determinations of the µ B -dependence of the QCD crossover transition line (see text). Right:
Net baryon-number density onthe lines of constant physics for three values of the energy density at µ B = 0. Other thermodynamic parameters characterizingthese lines are summarized in Table II. In the interval around T c , i.e. T ∈ [145 MeV ,
165 MeV] we find,0 . ≤ κ P ≤ . , . ≤ κ (cid:15) ≤ . , . ≤ κ s ≤ . . (50)Apparently, at O ( µ B ), lines of constant pressure and constant energy or entropy densities agree quite well and theyalso agree, within currently large errors, with the curvature of the transition line in (2+1)-flavor QCD. The coefficientof the quartic correction for the contour lines turns out to be about two orders of magnitude smaller than the leadingorder coefficients. This, of course, reflects the small contribution of the NLO corrections to the µ B -dependent partof pressure and energy density. For all fourth order coefficients we find | κ f | ≤ . T c . For µ B /T ≤ κ f only leads to modifications of T f ( µ B ) that stays withinthe error band arising from the uncertainty in κ f .The resulting lines of constant physics in the T - µ B plane are shown in Fig. 15 (left) for three values of the temper-ature, T = 145 MeV, 155 MeV and 165 MeV. These correspond to constant energy densities (cid:15) = 0 . ,0 . and 0 . , which roughly correspond to the energy density of cold nuclear matter, ahard sphere gas of nucleons at dense packing and the interior of a nucleus, respectively. Values of other bulk thermo-dynamic observables characterizing these LCPs are summarized in Table II. The corresponding net baryon-numberdensities on these LCPs are shown in Fig. 15 (right). It is apparent from Fig. 15 (left) that LCPs for constant pressure,1 at µ B = 0 on LCP T [MeV] p/T (cid:15)/T s/T p [GeV / fm ] (cid:15) [GeV / fm ] s [fm − ]145 0.586(80) 3.52(47) 4.11(53) 0.0337(46) 0.203(27) 1.63(21)155 0.726(95) 4.61(55) 5.34(63) 0.0546(71) 0.346(41) 2.59(30)165 0.898(110) 5.76(59) 6.66(69) 0.0868(106) 0.556(57) 3.90(40)Table II. Pressure, energy density and entropy density, characterizing lines of constant physics which correspond to the condi-tions met for µ B = 0 at T = 145 MeV, 155 MeV and 165 MeV. Columns 2-4 give results in appropriate units of temperature,while columns 5-7 give the same results expressed in units of GeV and fm . energy or entropy density agree well with each other up to baryon chemical potentials µ B /T = 2, where the differencein temperature on different LCPs is at most 2 MeV. We also note that the temperature on a LCP varies by about7 MeV or, equivalently, 5% between ˆ µ B = 0 and ˆ µ B = 2. Thus on a line of constant pressure, the entropy in units of T changes by about 15%. I.e. constant P or constant s/T , which both have been suggested as phenomenologicaldescriptions for freeze-out conditions in heavy ion collisions, can not hold simultaneously, although a change of 15%of one of these observables may phenomenologically not be of much relevance. We also stress that at large valuesof ˆ µ B the comparison of experimental data with HRG model calculations, e.g. the use of single particle Boltzmanndistributions used to extract freeze-out temperatures and chemical potentials, becomes questionable. As shown inFig. 11 net baryon-number densities extracted from HRG and QCD calculations differ substantially at µ B /T (cid:39) µ B /T ≤ √ s NN ≥ . T c ( µ B ) = T c (0) (cid:32) − κ c (cid:18) µ B T c (0) (cid:19) (cid:33) (51)with κ c ranging from 0 . . . . T c (0) = 155 MeV. While a small curvature for thecrossover line would suggest that the crossover transition happens under more or less identical bulk thermodynamicconditions a large curvature obviously would indicate that the crossover transition happens already at significantlysmaller values of pressure and energy density as µ B /T increases. VII. RADIUS OF CONVERGENCE AND THE CRITICAL POINT
As discussed in the previous sections we generally find that the Taylor series for all basic thermodynamic quantitiesconverge well for values of baryon chemical potentials µ B ≤ T . Even in the low temperature regime the relativecontribution of higher order expansion coefficients are generally smaller than in corresponding HRG model calculations.This, of course, also has consequences for our current understanding of the location of a possible critical point in theQCD phase diagram.The results on the expansion coefficients of the Taylor series for e.g. the pressure can be cast into estimates forthe location of a possible critical point in the QCD phase diagram. In general the radius of convergence can beobtained from ratios of subsequent expansion coefficients in the Taylor series for the pressure. Equally well one mayuse one of the derivatives of the pressure series. As one has to rely on estimates of the radius of convergence thatgenerally are based on a rather short series, it may indeed be of advantage to use as a starting point the series forthe net baryon-number susceptibility [47], which diverges at the critical point, but still contains information from allexpansion coefficients of the pressure series. The radius of convergence of this series is identical to that of the pressure.Model calculations also suggest that the estimators obtained from the susceptibility series converge faster to the trueradius of convergence [48]. For µ Q = µ S = 0 the expansion coefficients of the Taylor series for the net baryon-number2 r n χ -- e s t i m a t o r f o r µ B c r i t / T T [MeV]
Fodor, Katz, 2004Datta et al., 2016D’Elia et al., 2016, r χ this work: lower bound for r χ estimator r χ r χ ,HRG disfavored region for thelocation of a critical point r χ ,HRG Figure 16. Estimators for the radius of convergence of the Taylor series for net baryon-number fluctuations, χ B ( T, µ B ), in thecase of vanishing electric charge and strangeness chemical potentials obtained on lattices with temporal extent N τ = 8. Shownare lower bounds for the estimator r χ obtained in this work (squares) and results for this estimator obtained from calculationswith an imaginary chemical potential (triangles) [15]. Also shown are estimates for the location of the critical point obtainedfrom calculations with unimproved staggered fermions using a reweighting technique [50] and Taylor expansions [51]. In bothcases results have been rescaled using T c = 154 MeV. susceptibility are again simply related to that of the pressure, χ B ( T, µ B ) = ∞ (cid:88) n =0 n )! χ B n +2 ˆ µ nB . (52)From this one obtains estimators for the radius of convergence of the pressure and susceptibility series, r P n = (cid:12)(cid:12)(cid:12)(cid:12) (2 n + 2)(2 n + 1) χ B n χ B n +2 (cid:12)(cid:12)(cid:12)(cid:12) / , r χ n = (cid:12)(cid:12)(cid:12)(cid:12) n (2 n − χ B n χ B n +2 (cid:12)(cid:12)(cid:12)(cid:12) / . (53)Both estimators converge to the true radius of convergence in the limit n → ∞ . In order for this to correspond to asingularity at real values of ˆ µ B , all expansion coefficients should asymptotically stay positive.Obviously, the estimators r P n and r χ n are proportional to each other, r P n = (cid:112) (2 n + 2)(2 n + 1) / [2 n (2 n − r χ n .The difference between these to estimators may be taken as a systematic error for any estimate of the radius ofconvergence obtained from a truncated Taylor series. In the hadron resonance gas limit one finds for estimatorsinvolving sixth order cumulants, r P = 1 . r χ . In the following we restrict our discussion to an analysis of r χ n , whichat finite n leads to the smaller estimator for the radius of convergence. This seems to be appropriate in the presentsituation where we only can construct two independent estimators from ratios of three distinct susceptibilities. Wethus may hope to identify regions in the QCD phase diagram at small values of ˆ µ B which are unlikely locations for apossible critical point.An immediate consequence of the definitions given in Eq. 53 is that the ratios of generalized susceptibilities needto grow asymptotically like | χ Bn +2 /χ Bn | ∼ n in order to arrive in the limit n → ∞ at a finite value for the radiusof convergence. At least for large values of n one thus needs to find large deviations from the hadron resonancegas results | χ Bn +2 /χ Bn | HRG = 1. As is obvious from the results presented in the previous sections, in particularfrom Fig. 3, the analysis of up to sixth order Taylor expansion coefficients does not provide any hints for such largedeviations. The ratio χ B /χ B turns out to be less than unity in the entire temperature range explored so far, i.e.for T ≥
135 MeV or
T /T c > . T ∼
155 MeV, the sixth order expansioncoefficients also are consistent with HRG model results. They still have large errors. However, using the upper valueof the error for χ B /χ B provides a lower limit for the value of the estimator r χ . For temperatures in the interval135 MeV ≤ T ≤
155 MeV (or equivalently 0 . ≤ T /T c ≤
1) we currently obtain a lower limit on r χ from theestimate χ B /χ B (cid:39) χ B /χ B <
3. This converts into the bound r χ ≥
2, which is consistent with our observation thatthe Taylor series of all thermodynamic observables discussed in the previous sections is well behaved up to µ B = 2 T .A more detailed analysis, using the current errors on χ B /χ B at five temperature values below and in the crossoverregion of the transition at µ B = 0, is shown in Fig. 16. This shows that the bound arising from r χ is actually morestringent at temperatures closer to T c , where χ starts to become small and eventually tends to become negative.3These findings are consistent with recent results for susceptibility ratios obtained from calculations with an imag-inary chemical potential [15]. Also in that case all susceptibility ratios are consistent with HRG model results. Atpresent one thus cannot rule out that the radius of convergence may actually be infinite. Results for r χ obtained inRef. [15] lead to even larger estimators for the radius of convergence than our current lower bound. This is also shownin Fig. 16.The observations and conclusions discussed above are in contrast to estimates for the location of a critical pointobtained from a calculation based on a reweighting technique [50] as well as from Taylor series expansion in 2-flavorQCD [49, 51]. Both these calculations have been performed with unimproved staggered fermion discretization schemesand thus may suffer from large cut-off effects. Moreover, the latter calculation also suffers from large statistical errorson higher order susceptibilities. Results from Ref. [50] and Ref. [51] are also shown in Fig. 16.We thus conclude from our current analysis that a critical point at chemical potentials smaller than µ B = 2 T isstrongly disfavored in the temperature range 135 MeV ≤ T ≤
155 MeV and its location at higher values of temperatureseems to be ruled out. Our results suggest that the radius of convergence in that temperature interval will turn outto be significantly larger than the current bound once the statistics on 6 th order cumulants gets improved and higherorder cumulants become available. VIII. CONCLUSIONS
We have presented results on the equation of state of strong-interaction matter obtained from a sixth order Taylor-expansion of the pressure of (2+1)-flavor QCD with physical light and strange quark masses. We discussed expansionsat vanishing strangeness chemical potential µ S = 0 as well as for strangeness neutral systems n S = 0. We havediscussed in detail the latter case for a fixed electric charge to net baryon-number ratio, n Q /n B = 0 .
4, which isappropriate for situations met in heavy ion collisions. The results, however, can easily be extended to arbitrary ratiosof n Q /n B . We find that the dependence of basic thermodynamic observables on n Q /n B is small for 0 ≤ n Q /n B ≤ / T ∈ [130 MeV ,
330 MeV].We presented results for lines of constant pressure, energy and entropy density in the T - µ B plane and showed thatcorrections of O (ˆ µ B ) are negligible for ˆ µ B <
2. For all three observables the curvature term at O (ˆ µ B ) is smaller than κ max = 0 . κ max .The Taylor series for pressure and net baryon-number density as well as energy density and entropy density de-termined for µ S = 0 as well as n S = 0 have expansion coefficients that are close to HRG model results at lowtemperature. In general ratios of subsequent expansion coefficients approach the corresponding HRG model valuesfrom below when lowering the temperature. As a consequence, in the entire temperature range explored so far, theexpansions are ”better behaved” than the HRG model series, which have an infinite radius of convergence. Assumingthat the current results obtained with expansion coefficients up to 6 th order are indicative for the behavior of higherorder expansion coefficients and taking into account the current errors on 6 th order expansion coefficients we con-cluded that at temperatures T >
135 MeV the presence of a critical point in the QCD phase diagram for µ B ≤ T isunlikely. ACKNOWLEDGMENTS
This work was supported in part through Contract No. DE-SC001270 with the U.S. Department of Energy, throughthe Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy,Office of Science, Advanced Scientific Computing Research and Nuclear Physics, the DOE funded BEST topicalcollaboration, the NERSC Exascale Application Program (NESAP), the grant 05P12PBCTA of the German Bun-desministerium f¨ur Bildung und Forschung, the grant 56268409 of the German Academic Exchange Service (DAAD),grant 283286 of the European Union, the National Natural Science Foundation of China under grant numbers 11535012and 11521064 and the Early Career Research Award of the Science and Engineering Research Board of the Govern-ment of India. Numerical calculations have been made possible through an INCITE grant of USQCD, ALCC grantsin 2015 and 2016, and PRACE grants at CINECA, Italy, and the John von Neumann-Institute for Computing (NIC)in Germany. These grants provided access to resources on Titan at ORNL, BlueGene/Q at ALCF and NIC, CoriI and II at NERSC and Marconi at CINECA. Additional numerical calculations have been performed on USQCD4GPU and KNL clusters at JLab and Fermilab, as well as GPU clusters at Bielefeld University, Paderborn University,and Indiana University. We furthermore acknowledge the support of NVIDIA through the CUDA Research Centerat Bielefeld University.5
Appendix A: Details on simulation parameters and data sets
Our main data sets have been generated on lattices of size N σ × N τ , with N σ /N τ = 4 and N τ = 6 , m l /m s = 1 /
20 and 1 /
27, respectively.The simulation parameters are summarized in Table III and Table IV. N τ = 6 N τ = 8 β m l T[MeV] β m l T[MeV] m l /m s = 1 /
20 on lattices of size N σ N τ with N τ = 6 , N σ = 4 N τ . Columns 4 and 8 give the number of gauge field configurations, separated by 10 RHMC steps,that contributed to the analysis of up to sixth order generalized susceptibilities χ BQSijk . N τ = 6 N τ = 8 N τ = 12 β m l T[MeV] β m l T[MeV] β m l T[MeV] m l /m s = 1 /
27 and including results for N τ = 12. Appendix B: Constraints on chemical potential for strangeness neutral systems with fixed electric charge tobaryon-number ratio
We are interested in expansion coefficients for strangeness neutral systems in which the net electric-charge isproportional to the net baryon-number.
I.e. we introduce the constraint given in Eq. 31. These constraints can befulfilled order by order in the Taylor expansion of the number densities by choosing the expansion coefficients of theseries for ˆ µ Q and ˆ µ S , given in Eq. 14, appropriately, i.e. the coefficients s n and q n can be determined order by order.We start with the Taylor series for the number densities introduced in Eq. 16 and define the expansion coefficients as N Bn = s n χ BS + q n χ BQ + m Bn (B1) N Qn = s n χ QS + q n χ Q + m Qn (B2) N Sn = s n χ S + q n χ QS + m Sn (B3)6for n = 1 , ,
5. At each order in the expansion we then have to solve a set of two linear equations, which always havethe same structure. We find as solutions s n = − q n χ QS + m Sn χ S , (B4)and q n = − m Bn rχ S + m Qn χ S + m Sn ( rχ BS − χ QS )( χ QS ) − rχ BS χ QS + rχ S χ BQ − χ S χ Q . (B5)At leading order one finds for the terms m X , m B = χ B , m Q = χ BQ , m S = χ BS , (B6)and the contributions to the next-to-leading order expansion terms, m X , are given by m B = 16 (cid:0) q s χ BQS + 3 q s χ BQS + 6 q s χ BQS + q χ BQ + 3 q χ BQ +3 q χ BQ + s χ BS + 3 s χ BS + 3 s χ BS + χ B (cid:1) m Q = 16 (cid:0) q s χ QS + 3 q s χ QS + 6 q s χ BQS + q χ Q + 3 q χ BQ +3 q χ BQ + s χ QS + 3 s χ BQS + 3 s χ BQS + χ BQ (cid:1) m S = 16 (cid:0) q s χ QS + 3 q s χ QS + 6 q s χ BQS + q χ QS + 3 q χ BQS +3 q χ BQS + s χ S + 3 s χ BS + 3 s χ BS + χ BS (cid:1) (B7)Finally the contributions to the next-to-next-to-leading order expansion terms, m X , are given by m B = 1120 (cid:0) q s χ BQS + 10 q s χ BQS + 20 q s χ BQS + 60 q s χ BQS + 10 q s χ BQS + 30 q s χ BQS + 30 q s χ BQS +120 q s s χ BQS + 5 q s χ BQS + 120 q q s χ BQS + 120 q s χ BQS + 20 q s χ BQS + 30 q s χ BQS +20 q s χ BQS + 60 q s χ BQS + 120 q s χ BQS + q χ BQ + 5 q χ BQ + 10 q χ BQ + 60 q q χ BQ +10 q χ BQ + 120 q q χ BQ + 5 q χ BQ + 60 q χ BQ + 60 s s χ BS + s χ BS + 120 s s χ BS +5 s χ BS + 60 s χ BS + 10 s χ BS + 10 s χ BS + 5 s χ BS + χ B (cid:1) m Q = 1120 (cid:0) q s χ QS + 10 q s χ QS + 20 q s χ BQS + 60 q s χ QS + 10 q s χ QS + 30 q s χ BQS + 30 q s χ BQS +120 q s s χ QS + 5 q s χ QS + 120 q q s χ QS + 120 q s χ BQS + 20 q s χ BQS + 30 q s χ BQS +20 q s χ BQS + 60 q s χ QS + 120 q s χ BQS + q χ Q + 5 q χ BQ + 10 q χ BQ + 60 q q χ Q +10 q χ BQ + 120 q q χ BQ + 5 q χ BQ + 60 q χ BQ + 60 s s χ QS + s χ QS + 120 s s χ BQS +5 s χ BQS + 60 s χ BQS + 10 s χ BQS + 10 s χ BQS + 5 s χ BQS + χ BQ (cid:1) m S = 1120 (cid:0) q s χ QS + 10 q s χ QS + 20 q s χ BQS + 60 q s χ QS + 10 q s χ QS + 30 q s χ BQS + 30 q s χ BQS +120 q s s χ QS + 5 q s χ QS + 120 q q s χ QS + 120 q s χ BQS + 20 q s χ BQS + 30 q s χ BQS +20 q s χ BQS + 60 q s χ QS + 120 q s χ BQS + q χ QS + 5 q χ BQS + 10 q χ BQS + 60 q q χ QS +10 q χ BQS + 120 q q χ BQS + 5 q χ BQS + 60 q χ BQS + 60 s s χ S + s χ S + 120 s s χ BS +5 s χ BS + 60 s χ BS + 10 s χ BS + 10 s χ BS + 5 s χ BS + χ BS (cid:1) (B8)In (2+1)-flavor QCD calculations the light ( u, d ) quark masses are taken to be degenerate. A consequence of thisdegeneracy is that not all generalized susceptibilities χ BQSijk that enter the above expressions are independent. In a7 s T[MeV]
PDG-HRGQM-HRGcont. est.N τ =161286
20 (open)m s /m l =27 (filled) free quark gasn S =0, n Q /n B =0.4 - q T[MeV]
PDG-HRGQM-HRGcont. est.N τ =161286
20 (open)m s /m l =27 (filled) free quark gasn S =0, n Q /n B =0.4 s / s T[MeV]
PDG-HRGQM-HRGcont. est.N τ =86
20 (open)m s /m l =27 (filled)free quark gasn S =0, n Q /n B =0.4 -0.02 0 0.02 0.04 0.06 0.08 0.1 140 160 180 200 220 240 260 280 q / q T[MeV]
PDG-HRGQM-HRGcont. est.N τ =86 -0.02 0 0.02 0.04 0.06 0.08 0.1 140 160 180 200 220 240 260 280
20 (open)m s /m l =27 (filled)free quark gasn S =0, n Q /n B =0.4 Figure 17. The LO Taylor expansion coefficients s (top, left) and q (top, right) of the expansions of ˆ µ S and ˆ µ Q w.r.t. ˆ µ B . Thebottom set of figures show the ratios of NLO and LO expansion coefficients. The broad bands give the continuum extrapolatedresults. The curves inside these bands show results obtained with the interpolating curves introduced in Eq. 38. Also shown arethe PDG-HRG and QM-HRG results (see text). The solid black lines labeled ‘free quark gas’ denote the T → ∞ non-interactingmassless quark gas result. given order n ≡ l ≡ i + j + k this results in a set of relations among the expansion coefficients. In general, at order n = 2 l , there are l ( l + 1) constraints, i.e. for l = 1 this gives rise to two relations, [27]0 = χ B − χ BQ + χ BS χ S − χ QS + χ BS , (B9)for l = 2 there are six constraints,0 = χ B − χ BQ + χ BS χ S − χ QS + χ BS χ BS + χ BS − χ BQS χ BS + χ BS − χ BQS χ B − χ BQ + 12 χ BQ − χ BQ + 3 χ BS + 3 χ BS + χ BS − χ BQS + 12 χ BQS − χ BQS χ S + χ BS + 3 χ BS + 3 χ BS − χ QS + 12 χ QS − χ QS − χ BQS + 12 χ BQS − χ BQS (B10)8and for l = 3 there are twelve constraints,0 = χ B − χ BQ + χ BS , χ BS − χ QS + χ S , χ BS − χ BQS + χ BS , χ BS − χ BQS + χ BS , χ BS − χ BQS + χ BS , χ BS − χ BQS + χ BS , χ B − χ BQ + 12 χ BQ − χ BQ + 3 χ BS − χ BQS + 12 χ BQS + 3 χ BS − χ BQS + χ BS , χ BS − χ BQS + 12 χ BQS − χ QS + 3 χ BS − χ BQS + 12 χ QS + 3 χ BS − χ QS + χ S , χ BS − χ BQS + 12 χ BQS − χ BQS + 3 χ BS − χ BQS + 12 χ BQS + 3 χ BS − χ BQS + χ BS , χ BS − χ BQS + 12 χ BQS − χ BQS + 3 χ BS − χ BQS + 12 χ BQS + 3 χ BS − χ BQS + χ BS , χ B − χ BQ + 40 χ BQ − χ BQ + 80 χ BQ − χ BQ + 5 χ BS − χ BQS + 120 χ BQS − χ BQS + 80 χ BQS + 10 χ BS − χ BQS + 120 χ BQS − χ BQS + 10 χ BS − χ BQS +40 χ BQS + 5 χ BS − χ BQS + χ BS , χ BS − χ BQS + 40 χ BQS − χ BQS + 80 χ BQS − χ QS + 5 χ BS − χ BQS + 120 χ BQS − χ BQS + 80 χ QS + 10 χ BS − χ BQS + 120 χ BQS − χ QS + 10 χ BS − χ BQS +40 χ QS + 5 χ BS − χ QS + χ S . (B11)Using these constraints it is tedious, but straightforward, to show that in the isospin symmetric case, r = 1 /
2, indeedall expansion coefficients for the electric charge chemical potential vanish, i.e. ˆ µ Q = 0 to all orders in µ B .We show results for the LO expansion coefficients s and q and the ratios of the NLO and LO expansion coefficients, s /s and q /q in Fig. 17. As can be seen the NLO coefficients are already negligible for T > ∼
170 MeV. The absolutevalue of the NNLO expansion coefficients s and q never is larger than 1% of the corresponding LO coefficients.In Fig. 17, we also show results from hadron resonance gas (HRG) model calculations. The black curves are thepredictions of the usual HRG model which consists of all the resonances listed in the Particle Data Group Tables up to2.5 GeV (PDG-HRG). The PDG-HRG results for s are substantially smaller than the continuum extrapolated latticeQCD results. It has been argued in [6] that this can be caused by contributions from additional, experimentally notyet observed, strange hadron resonances which are predicted in quark model calculations. A HRG model calculationbased on such an extended resonance spectrum (QM-HRG) is also shown in Fig. 17. At finite values of the latticecut-off we observe significant differences between lattice QCD calculations and both versions of the HRG models. Thisis in particular the case for the expansion coefficients of the electric charge chemical potentials. One thus may wonderwhether these deviations can be understood in terms of taste violations in the staggered fermion formulation whichresult in a modification of the resonance spectrum and affect most strongly the light pseudo-scalar (pion) sector. Appendix C: The coefficient κ f of lines of constant physics at O ( µ B ) We will present here results for the expansion coefficient κ f of lines of constant physics defined in Eq. (43), κ f = ∂ f ( T,µ B ) ∂T (cid:12)(cid:12)(cid:12) ( T , ( κ f ) T − ∂∂T ∂ f ( T,µ B ) ∂µ B (cid:12)(cid:12)(cid:12) ( T , κ f T + ∂ f ( T,µ B ) ∂µ B (cid:12)(cid:12)(cid:12) ( T , ∂f ( T,µ B ) ∂T (cid:12)(cid:12)(cid:12) ( T ,
0) 1 T = T ∂ f ( T ) ∂T (cid:12)(cid:12)(cid:12) ( T , ( κ f ) − (cid:16) T ∂f ( T ) ∂T (cid:12)(cid:12)(cid:12) ( T , − f ( T ) (cid:19) κ f + f ( T ) T ∂f ( T ) ∂T (cid:12)(cid:12)(cid:12) ( T , (C1)The coefficients f k are defined by f ( T, µ B ) = ∞ (cid:88) k =0 f k ˆ µ kB . (C2)9In particular, we will give explicit expressions for the case of constant pressure ( f ≡ P ), constant energy density( f ≡ (cid:15) ) and constant entropy density( f ≡ s ). For the pressure we had the earlier expression (Eq. (15)) P ( T, µ B ) − P ( T, T = ∞ (cid:88) n =1 P n ˆ µ nB . (C3)Comparing Eqs. (C3) and (C2) we have, f = P ( T, ≡ T P , f = T P and f = T P . Thus, ∂f ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B = ∂P T ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B = T ( T P (cid:48) + 4 P ) ≡ s, (C4a) ∂ f ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B = ∂ P T ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B = T (cid:0) T P (cid:48)(cid:48) + 8 T P (cid:48) + 12 P (cid:1) ≡ C V T . (C4b)Here s and C V are the entropy density and specific heat per unit volume at vanishing chemical potential. Similarly, ∂f ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B = ∂P T ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B = T ( T P (cid:48) + 4 P ) . (C5)Putting everything together we get, for the pressure: κ P = 1 T P (cid:48) + 4 P (cid:20) P − κ P ( T P (cid:48) + 2 P ) + 12 (cid:0) κ P (cid:1) (cid:0) T P (cid:48)(cid:48) + 8 T P (cid:48) + 12 P (cid:1)(cid:21) = T s (cid:20) P ( T ) − κ P σ ( T ) + 12 (cid:0) κ P (cid:1) C V T (cid:21) , (C6)where σ denotes the O (ˆ µ B ) expansion coefficient of the entropy density as introduced in Eq. 24.Next we consider κ (cid:15) . Since the energy density is also of dimension four, we only need to replace P n with (cid:15) n inthe first line of Eq. (C6). With this we obtain, κ (cid:15) = 1 T (cid:15) (cid:48) + 4 (cid:15) (cid:20) (cid:15) − κ (cid:15) ( T (cid:15) (cid:48) + 2 (cid:15) ) + 12 ( κ (cid:15) ) (cid:0) T (cid:15) (cid:48)(cid:48) + 8 T (cid:15) (cid:48) + 12 (cid:15) (cid:1)(cid:21) . (C7)Since C V ≡ ( ∂(cid:15) /∂T ) µ B , the above may be written as κ (cid:15) = T C V (cid:20) (cid:15) − κ (cid:15) ( T (cid:15) (cid:48) + 2 (cid:15) ) + 12 ( κ (cid:15) ) T ∂C V ∂T (cid:21) . (C8)Finally we consider κ s . Since the entropy density is of dimension three, Eqs. (C4) become ∂ ( sT ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B = T ( T s (cid:48) + 3 s ) , ∂ ( sT ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) µ B = T (cid:0) T s (cid:48)(cid:48) + 6 T s (cid:48) + 6 s (cid:1) , (C9)and therefore κ σ = 1 T s (cid:48) + 3 s (cid:20) σ − κ σ ( T σ (cid:48) + σ ) + 12 ( κ σ ) (cid:0) T s (cid:48)(cid:48) + 6 T s (cid:48) + 6 s (cid:1)(cid:21) . (C10)To zeroth order, the specific heat is also given by C V = ( ∂ ( T s ) /∂T ) µ B . Thus, κ σ = T C V (cid:20) σ − κ σ ( T σ (cid:48) + σ ) + 12 ( κ σ ) T ∂C V ∂T (cid:21) . (C11) [1] J. Engels, F. Karsch, H. Satz and I. Montvay, Phys. Lett. B , 89 (1981).[2] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabo, Phys. Lett. B
99 (2014) [arXiv:1309.5258[hep-lat]]. [3] A. Bazavov et al. [HotQCD Collaboration], Phys. Rev. D , 094503 (2014) [arXiv:1407.6387 [hep-lat]].[4] A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. T. Ding, S. Gottlieb, R. Gupta and P. Hegde et al. , Phys. Rev. D , 252002 (2010) [arXiv:1008.1747 [hep-ph]].[6] A. Bazavov, H.-T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa and S. Mukherjee et al. , Phys.Rev. Lett. , 072001 (2014) [arXiv:1404.6511 [hep-lat]].[7] R. V. Gavai and S. Gupta, Phys. Rev. D , 034506 (2003) [hep-lat/0303013].[10] C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and C. Schmidt, Phys. Rev. D , 014507(2003) [hep-lat/0305007].[11] C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and K. Redlich, Phys. Rev. D ,054508 (2005) [hep-lat/0501030].[12] S. Ejiri, F. Karsch, E. Laermann and C. Schmidt, Phys. Rev. D , 054506 (2006) [hep-lat/0512040].[13] S. Borsanyi, G. Endrodi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti and K. K. Szabo, JHEP (2012) 053 [arXiv:1204.6710[hep-lat]].[14] J. Gunther, R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz, A. Pasztor and C. Ratti, arXiv:1607.02493 [hep-lat].[15] M. D’Elia, G. Gagliardi and F. Sanfilippo, arXiv:1611.08285 [hep-lat].[16] M. Asakawa, K. Yazaki, Nucl. Phys. A , 668 (1989).[17] A. M. Halasz et al., Phys. Rev. D , 096007 (1998).[18] A. Vuorinen, Phys. Rev. D , 054017 (2003) doi:10.1103/PhysRevD.68.054017 [hep-ph/0305183].[19] E. Follana et al. [HPQCD and UKQCD Collaborations], Phys. Rev. D et al. [HotQCD Collaboration], Phys. Rev. D , 308 (1983).[22] R. V. Gavai and S. Sharma, Phys. Rev. D , 054508 (2012) doi:10.1103/PhysRevD.85.054508 [arXiv:1112.5428 [hep-lat]].[23] R. V. Gavai and S. Sharma, Phys. Lett. B , 8 (2015) doi:10.1016/j.physletb.2015.07.036 [arXiv:1406.0474 [hep-lat]].[24] B. Friman, F. Karsch, K. Redlich and V. Skokov, Eur. Phys. J. C , 1694 (2011) doi:10.1140/epjc/s10052-011-1694-2[arXiv:1103.3511 [hep-ph]].[25] C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro and F. Sanfilippo, Phys. Rev. D , no. 7, 074504 (2016)doi:10.1103/PhysRevD.93.074504 [arXiv:1602.01426 [hep-lat]].[26] A. Roberge and N. Weiss, Nucl. Phys. B , 734 (1986).[27] A. Bazavov, H. T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, S. Mukherjee and P. Petreczky et al. , Phys.Rev. Lett. , 192302 (2012) [arXiv:1208.1220 [hep-lat]].[28] R. Hagedorn and J. Rafelski, Phys. Lett. , 136 (1980).[29] V. V. Dixit, F. Karsch and H. Satz, Phys. Lett. , 412 (1981).[30] A. Andronic, P. Braun-Munzinger, J. Stachel and M. Winn, Phys. Lett. B , 80 (2012) [arXiv:1201.0693 [nucl-th]].[31] V. Vovchenko, M. I. Gorenstein and H. Stoecker, arXiv:1609.03975 [hep-ph].[32] J. Cleymans and K. Redlich, Phys. Rev. C , 054908 (1999).[33] J. Cleymans, H. Oeschler, K. Redlich and S. Wheaton, Phys. Rev. C , 034905 (2006).[34] B. Tomasik and U. A. Wiedemann, Phys. Rev. C , 034905 (2003) [nucl-th/0207074].[35] J. Cleymans et al. [NA49 Collaboration], Phys. Lett. B , 50 (2005) [hep-ph/0411187].[36] J. Rafelski and J. Letessier, J. Phys. G , 064017 (2009) [arXiv:0902.0063 [hep-ph]].[37] M. Petran and J. Rafelski, Phys. Rev. C , no. 2, 021901 (2013).[38] J. Rafelski and M. Petran, Phys. Part. Nucl. , no. 5, 748 (2015).[39] S. Das [STAR Collaboration], EPJ Web Conf. , 08007 (2015) [arXiv:1412.0499 [nucl-ex]].[40] M. Floris, Nucl. Phys. A , 103 (2014) [arXiv:1408.6403 [nucl-ex]].[41] F. Becattini, J. Steinheimer, R. Stock and M. Bleicher, Phys. Lett. B , 241 (2017) [arXiv:1605.09694 [nucl-th]].[42] O. Kaczmarek, F. Karsch, E. Laermann, C. Miao, S. Mukherjee, P. Petreczky, C. Schmidt and W. Soeldner et al. , Phys.Rev. D , 014504 (2011) [arXiv:1011.3130 [hep-lat]].[43] G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, JHEP , 001 (2011) [arXiv:1102.1356 [hep-lat]].[44] C. Bonati, M. D’Elia, M. Mariti, M. Mesiti, F. Negro and F. Sanfilippo, Phys. Rev. D , no. 5, 054503 (2015)doi:10.1103/PhysRevD.92.054503 [arXiv:1507.03571 [hep-lat]].[45] R. Bellwied, S. Borsanyi, Z. Fodor, J. G ˜A nther, S. D. Katz, C. Ratti and K. K. Szabo, Phys. Lett. B , 559 (2015)doi:10.1016/j.physletb.2015.11.011 [arXiv:1507.07510 [hep-lat]].[46] P. Cea, L. Cosmai and A. Papa, Phys. Rev. D , no. 1, 014507 (2016) doi:10.1103/PhysRevD.93.014507 [arXiv:1508.07599[hep-lat]].[47] R. V. Gavai and S. Gupta, Phys. Rev. D , 114014 (2005) doi:10.1103/PhysRevD.71.114014 [hep-lat/0412035].[48] F. Karsch, B. J. Schaefer, M. Wagner and J. Wambach, PoS LATTICE , 219 (2011) [arXiv:1110.6038 [hep-lat]].[49] S. Datta, R. V. Gavai and S. Gupta, PoS LATTICE , 202 (2014).[50] Z. Fodor and S. D. Katz, JHEP0404