Lattice QCD-based equations of state at vanishing net-baryon density

We present realistic equations of state for QCD matter at vanishing net-baryon density which embed recent lattice QCD results at high temperatures combined with a hadron resonance gas model in the low-temperature, confined phase. In the latter, we allow an implementation of partial chemical equilibrium, in which particle ratios are fixed at the chemical freeze-out, so that a description closer to the experimental situation is possible. Given the present uncertainty in the determination of the chemical freeze-out temperature from first-principle lattice QCD calculations, we consider different values within the expected range. The corresponding equations of state can be applied in the hydrodynamic modeling of relativistic heavy-ion collisions at the LHC and at the highest RHIC beam energies. Suitable parametrizations of our results as functions of the energy density are also provided.


I. INTRODUCTION
In the relativistic heavy-ion collisions at RHIC (Relativistic Heavy-Ion Collider) and LHC (Large Hadron Collider), a hot deconfined state of strongly interacting matter is transiently created, the Quark-Gluon Plasma (QGP). This form of QCD matter is believed to have existed in the very first moments of our universe. As the produced hot and dense system cools down during its expansion, matter undergoes a transition from the QGP phase into a state dominated by color-confined, massive hadronic degrees of freedom. The nature of this phase transformation has been determined at vanishing baryo-chemical potential by first-principle lattice QCD simulations: it is an analytic crossover, taking place over a broad region of temperatures T [1]. The value of the (pseudo-) critical temperature T c associated with this confinement transition depends to some extent on the considered order-parameter. For example, the Wuppertal-Budapest (WB) and hotQCD collaborations found comparable values for chiral symmetry restoration: T c =(155 ± 6) MeV in [2] and T c =(154 ± 9) MeV in [3], respectively.
The collective flow dynamics of the bulk of matter created in heavy ion collisions can be successfully modeled by means of relativistic hydrodynamics (cf. e.g. the reviews in [4,5]), starting from a stage immediately after thermalization until the kinetic freeze-out of final state hadrons. Assuming local thermal equilibrium, the conservation equations for energy, momentum and for the additionally conserved charges (net-baryon number N B , net-electric charge N Q and net-strangeness N S ) drive the evolution of the system. An essential ingredient for this modeling is the equation of state (EoS), which provides locally a relation between energy density ǫ, pressure p * Electronic address: mbluhm@to.infn.it and the densities n B , n Q and n S of the conserved charges. The parameter controlling the acceleration of the fluid's collective flow by pressure gradients is the speed of sound, c s = ∂p/∂ǫ.
A quantitative comparison of hydrodynamic simulations with the observed collective flow behavior revealed that the evolution of the system can be described by nearly ideal hydrodynamics, cf. e.g. [6][7][8][9][10][11][12]. In these studies, a uniquely small ratio of shear viscosity η to entropy density s of the hot matter was determined, cf. also the reviews in [13,14]. This led to our current understanding of the QGP as a strongly coupled, nearly perfect fluid [15][16][17]. Assuming the conservation of entropy, i.e. neglecting the viscous entropy production associated with such a small η/s [18], one needs to know the EoS only along adiabatic paths. In this work, we concentrate on the situation of a vanishing n B , i.e. we consider the path n B /s = 0. We note that in the thermal system created in a heavy-ion collision one always has n S = 0, while n Q (in the case of a partial stopping at the lowest center-of-mass energies) is related to n B .
A rigorous determination of the equation of state for n B = 0 in the non-perturbative regime of QCD can be achieved with lattice gauge theory simulations. These reach nowadays unprecedented levels of accuracy. A basic quantity for the EoS is the interaction measure I = ǫ−3p, which has been calculated in [19,20] and in [21][22][23]. The numerical results for I(T )/T 4 in [19,20] show significant differences from those in [22,23] in the transition region. We opt for utilizing in our work the most recent, continuum-extrapolated lattice QCD data from the WB-collaboration [23], corresponding to a system of 2+1 quark flavors with physical quark masses.
By combining a suitable parametrization of these lattice QCD results with a hadron resonance gas (HRG) model in chemical equilibrium, we construct a baseline QCD equation of state for n B = 0. Moreover, we consider also the case of a HRG in partial chemical equilibrium in order to properly account for the actual chem-ical composition in the hadronic phase. This is known to be of importance in order to reproduce not only the observed flow and p T -spectra, but also the correct particle ratios [24]. Because of the present uncertainty in its exact value [25], we study various values for the chemical freeze-out temperature T ch , below which the HRG is in partial chemical equilibrium. In this way, different QCD equations of state are obtained, which can be used in the hydrodynamic simulations of relativistic heavyion collisions for LHC and RHIC top beam energies at mid-rapidity when net-baryon density effects can be neglected. Having such QCD equations of state at hand will allow a more controlled determination of the QGP transport properties, as for example the shear (and bulk) viscosity coefficients. At smaller beam energies, effects of a non-vanishing net-baryon density become important. Corresponding QCD equations of state will be presented in a forthcoming publication.
The equation of state of QCD matter has been the subject of numerous studies in the literature. Among different other approaches, we mention combinations of the HRG model with an effective theory of QCD [26], with a phenomenological model for QCD thermodynamics [27] and with various parametrizations [28][29][30][31] of lattice QCD results. Developments in using a parametrization of lattice QCD results for finite n B were recently reported in [32,33]. Moreover, in [34] an EoS, describing both the QGP and the hadronic phase based on one effective model approach, was constructed and applied in finite-n B hydrodynamics studies (see also further developments in [35]).
The paper is organized as follows: in section II, we discuss briefly the employed lattice QCD results [23] and their combination with a HRG model in chemical equilibrium. Section III deals with the inclusion of partial chemical equilibrium in the description of the hadronic phase. In section IV, we discuss the obtained QCD equations of state and provide practical parametrizations of our results.

II. CONSTRUCTION OF A LATTICE QCD-BASED EOS
In [23], continuum-extrapolated lattice gauge theory results of QCD thermodynamics for 2+1 quark flavors with physical mass parameters were presented. The corresponding results for the scaled interaction measure I(T )/T 4 and for the scaled pressure p(T )/T 4 are depicted in Fig. 1 panels (a) and (b), respectively. The other thermodynamic quantities are related to I(T ) and p(T ) via ǫ(T ) = 3p(T )+I(T ) and s(T ) = (ǫ(T )+p(T ))/T . A suitable parametrization of these lattice QCD results yields for s(T ) extrapolated to T = 800 MeV a value of about 82.5% of the Stefan-Boltzmann limit of a non-interacting gas of 3 massless quark flavors.
The QCD thermodynamics in the hadronic phase can be well accounted for by the HRG model describing hadronic matter in thermal and chemical equilibrium, cf. e.g. [36,37]. The pressure of the model in the thermodynamic limit is given by where the sum is taken over all hadronic (including resonances) states k (baryons and anti-baryons being summed independently) included in the model. In Eq. (1), d k and m k denote the degeneracy factor and the mass, and µ k is the chemical potential of the hadronspecies k. In chemical equilibrium, the latter reads where B k , Q k and S k are the respective quantum numbers of baryon charge, electric charge and strangeness, while µ B , µ Q and µ S denote the chemical potentials associated with n B , n Q and n S .
Other thermodynamic quantities follow from standard relations, e.g. s = (∂p/∂T ) µ k . The particle number density of species k, n k = (∂p/∂µ k ) T , is given by the momentum-integral (2) and the net-baryon density follows from n B = k B k n k . Since we consider n B = 0, all µ k are set to zero in the chemical equilibrium case.
In this work, we employ a HRG model containing states up to a mass of 2 GeV as, for example, listed in the edition [38] of the Particle Data Book. Such a list is also included in the EoS-package provided along with the work in [31]. As evident from Fig. 1, this choice is enough to describe the available lattice QCD data fairly well for temperatures below 175 MeV, where HRG and lattice QCD results overlap. In fact, the relative deviation of the HRG model from the lattice QCD data [23] in this overlap region, taking the data error-bars into account, is at most 9% in I(T )/T 4 and 5% in p(T )/T 4 .
Given the reasonable agreement between lattice QCD data and the HRG model, we construct an equation of state, which may serve as a baseline: we utilize the lattice QCD results from [23] at high T and change the prescription to a HRG model at low T around a switching temperature of 172 MeV. Generically, such an approach can introduce discontinuities in the thermodynamic quantities. We improve this situation by employing a straightforward interpolation procedure in the interval 165 MeV ≤ T ≤ 180 MeV, which ensures that the pressure and its first and second derivatives with respect to T are continuous. In this way, the speed of sound remains a smooth function for all temperatures.

III. HADRON RESONANCE GAS IN PARTIAL CHEMICAL EQUILIBRIUM
In heavy-ion collisions, the time scales for inelastic particle number changing processes, which are responsible for the chemical equilibration of the hadronic matter, are typically much larger than the lifetime of the hadronic stage [39]. Thus, it is more reasonable to assume that the hadronic phase is not in complete chemical equilibrium. This was first discussed in [40] and then considered in numerous works, cf. e.g. [24,[41][42][43][44]: according to this idea, hadronic matter is formed at the hadronization temperature T c in chemical equilibrium. However, for temperatures below the chemical freeze-out temperature T ch , where T ch ≤ T c , the inelastic processes become suppressed, while the elastic interactions mediated by frequent strong resonance formations and decays (e.g. ππ → ρ → ππ, Kπ → K * → Kπ, pπ → ∆ → pπ etc) continue to occur. Consequently, the experimentally observed ratios of particle multiplicities of those species i, which are stable against strong decays within the lifetime of the system, are fixed at T ch . This is to say that for T < T ch the corresponding effective particle num-bersN i = N i + r d r→i N r are frozen. Here, N i denotes the actual particle number of the stable hadron i, N r the actual particle number of resonance r and d r→i gives the average number of hadrons i produced in the decay of resonance r. For example, the conserved quantity in the process ππ → ρ → ππ is the effective pion number N π = N π + 2N ρ . The above sum has to be taken over all the states (resonances) that decay into hadron i within the lifetime of the hadronic stage. As their effective number is fixed at T ch , but T decreases during the expansion of matter, each stable particle species i acquires an effective, T -dependent chemical potential µ i (T ). The chemical potentials of the resonances, instead, can be written as a combination µ r = i d r→i µ i of the effective chemi-cal potentials. The hadronic phase is, thus, in a state of partial chemical equilibrium below T ch .
The freeze-out of the chemical composition of the system at T ch implies, in addition to the conservation of energy, momentum and of the charges N B , N Q and N S , also the conservation of the effective numberN i of each stable particle species i below T ch . This makes the EoS a highly-involved relation between p, ǫ and all charge densities. Assuming the conservation of entropy, the ratio between the effective particle number density and the entropy densityn i /s is fixed at T ch . This provides a practical tool to conserve all theN i and to determine all the µ i (T ) for T < T ch from the conditions which imply that eachn i depends, in general, on all the effective chemical potentials µ i ′ (T ) (including µ i (T )). The knowledge of all the µ i (T ) is, apart from knowing the EoS, necessary for determining the final state hadron abundances. We note that the above conditions entail also that the particle ratios of stable hadrons are fixed at T ch :n i1 (T, In this work, we consider in line with [31] as stable particle species the mesons π 0 , π + , π − , K + , K − , K 0 , K 0 and η and the baryons p, n, Λ 0 , Σ + , Σ 0 , Σ − , Ξ 0 , Ξ − and Ω − as well as their respective anti-baryons, i.e. in total 26 different states. Correspondingly, we consider different isospin states individually. In general, this becomes important only when considering non-vanishing net-densities n B , n Q and/or n S . In the n B = 0 case studied in this work, however, particles and their corresponding anti-particles develop the same effective chemical potentials. For the chemical freeze-out temperature, we consider different values, namely T ch /MeV=    145, 150, 155 and 160. These are within the range of the T c -values determined in lattice QCD [2,3]. Moreover, we do not choose larger T ch -values, because otherwise contributions from resonance states heavier than 2 GeV cannot be neglected anymore in the analysis. In Fig. 2, we exhibit the temperature-dependence of the effective chemical potentials µ i (T ) of some exemplary particle species as determined from Eq. (3) for T ch = 150 MeV. As can be seen from Fig. 2, the µ i (T ) increase with decreasing T . The T -dependence of µ i (T ) for species i can be parametrized conveniently by the quadratic function Here, the parameters a i and b i depend on the value of T ch . Since for a complete EoS the knowledge of all µ i (T ) is required, we summarize the corresponding parametervalues in Tab. I. With these, µ i (T ) is obtained in units of GeV for T ch and T given in units of GeV. These parametrizations provide excellent fits for all µ i (T ) in the temperature range 70 MeV ≤ T ≤ T ch with a maximal χ 2 = 9 · 10 −6 . We note that overall, for the large T -range explored in a hydrodynamic simulation, cubic functions for the µ i (T ) yield more accurate descriptions of the numerical results obtained from Eq. (3) than the quadratic functions, in particular for small T . In the interesting interval 70 MeV ≤ T ≤ T ch , however, the quadratic ansatz Eq. (4) provides fits, which are comparable in accuracy with the cubic-fits for the baryons and anti-baryons, while they are even slightly better for the mesons.

IV. DISCUSSION AND CONCLUSIONS
Combining the lattice QCD data [23] with the HRG model in chemical equilibrium and then considering optionally partial chemical equilibrium in the hadronic phase with various T ch -values, we obtain different QCD equations of state. For the use in a hydrodynamic simulation, the EoS is usually given in the form p(ǫ, n B ) together with the results for the effective chemical potentials µ i required for determining the particle abundances. In Fig. 3, we show our results for the different equations of state p(ǫ) supplemented by T (ǫ) for n B = 0. We concentrate in Fig. 3  density regions, in which the confinement transition and the chemical freeze-out take place.
As it is evident from Fig. 3 panel (a), differences in p(ǫ) between chemical equilibrium (solid curve) and partial chemical equilibrium (dashed curves) in the hadronic phase are small. The ǫ-dependence of T , instead, is significantly influenced for ǫ < ǫ ch by the chemical freeze-out (cf. panel (b) in Fig. 3), where the value of ǫ ch depends on T ch .
Our results are collected in tabulated form and made available along with this publication [45]. Moreover, for convenience we also provide suitable parametrizations of these numerical results, similar to Ref. [18]. When assuming chemical equilibrium in the hadronic phase, the relevant thermodynamic quantities can be parametrized in the following way: and Here, Γ(s, x) = ∞ x t s−1 exp[−t] dt denotes the upper incomplete Γ-function. These ansatz-functions can provide excellent descriptions of our numerical results for proper choices of the entering parameters. We stress, that the parameters a i in p(ǫ) and s 4/3 (ǫ) in Eqs. (5) and (6) are not meant to be the same: we use the same symbols for practical purposes.
It turns out that, for an accurate description of the thermodynamic quantities, it becomes mandatory to split the parametrizations into different regions in ǫ and to fit the parameters for each ǫ-region individually. We define as the splitting-points ǫ 1 = 0.032084 GeV/fm 3 , ǫ 2 = 0.567420 GeV/fm 3 and ǫ 3 = 100 GeV/fm 3 . The latter point ǫ 3 is of relevance for the parametrization of s 4/3 (ǫ) in Eq. (6) and therefore influences T (ǫ), but plays no role for p(ǫ). The parameter-values entering p(ǫ) and s 4/3 (ǫ) in the different ǫ-regimes are summarized in Tab. II. With these, one obtains p in units of GeV/fm 3 , s in units of 1/fm 3 and T from Eq. (7) in units of GeV for ǫ in units of GeV/fm 3 . We make sure that thermodynamic consistency in all quantities is maintained at the splitting-points ǫ i up to a numerical accuracy of 1 · 10 −8 .
The squared speed of sound c 2 s (ǫ) can be determined from p(ǫ) given in Eq. (5) as By employing the parameter-values for p(ǫ) from Tab. II, we find an excellent agreement between Eq. (8) and the c 2 s -results obtained from a numerical differentiation of our p(ǫ)-results for all ǫ > ǫ 1 with an accuracy of 1 · 10 −5 or better. For ǫ < ǫ 1 , the agreement between Eq. (8) and the numerical results is still reasonable, although for ǫ < 0.0008 GeV/fm 3 the fit Eq. (8) (5) and (6)  For comparison, the symbols depict available equilibrium lattice QCD data from [22]. The dashed curves highlight c 2 s (T ) when instead partial chemical equilibrium is assumed in the hadronic phase. We consider T ch /MeV = 145, 150, 155 and 160 (from top to bottom, respectively).
approaches no unphysical result in the limit of ǫ → 0.
The temperature-dependence of c 2 s is shown in Fig. 4 (solid curve) and confronted with the available lattice QCD results of the WB-collaboration reported in [22]. In line with these lattice data, we find a rather large c 2 s in the confinement transition region. This indicates that our EoS is rather stiff compared to some previously considered equations of state, as e.g. in [18].
When including partial chemical equilibrium into the EoS, the parametrizations discussed above become modified only for ǫ < ǫ ch . The different values of ǫ ch , depending on the considered value for the chemical freeze-out temperature, are listed in the caption of Tab. III. For ǫ < ǫ ch , we modify our parametrizations to and Correspondingly, the squared speed of sound follows now from Eq. (9) as We stress that the parameters b i entering p(ǫ), s 4/3 (ǫ) and T (ǫ) in Eqs. (9) -(11) are also here not meant to be the same. By fitting the parametrizations in Eqs. (9) -(11) to our numerical results for ǫ < ǫ ch , we find quite accurate descriptions of all thermodynamic quantities, including c 2 s , for the parameter-values summarized in Tab. III. With these, p, s and T are obtained in units of GeV/fm 3 , 1/fm 3 and GeV, respectively, for ǫ given in units of GeV/fm 3 . Moreover, the above parametrizations satisfy the physical conditions T (ǫ) → 0, p(ǫ) → 0 and s(ǫ) → 0 for ǫ → 0 as well as 0 ≤ c 2 s ≤ 1/3 (indeed c 2 s (ǫ → 0) = b 3 b 4 ). The temperature-dependence of c 2 s as determined from a numerical differentiation is shown in Fig. 4 (dashed curves). One observes a discontinuity in c 2 s (T ) at T = T ch , which is characteristic for the chemical freeze-out.
In summary, we constructed QCD equations of state for vanishing net-baryon density based on recent continuum-extrapolated lattice QCD results in the physical quark mass limit at high T [23], continuously combined with a HRG model at low T . The latter was considered to be either in chemical equilibrium or in partial chemical equilibrium. In the partial chemical equilibrium case, we studied different values for the chemical freezeout temperature within its presently expected range [25].
Our results, being available in a tabulated form [45], can be directly applied in the hydrodynamic modeling of high-energy heavy-ion collisions at the LHC and at RHIC for top beam energies and at mid-rapidity. For convenience, we also provided suitable parametrizations of our   (10) and (11) for p(ǫ), s 4/3 (ǫ) and T (ǫ), respectively, providing suitable parametrizations of our results for ǫ < ǫ ch when partial chemical equilibrium is considered in the hadronic phase. For ǫ > ǫ ch , the thermodynamic quantities are given by Eqs. results, in particular, for the effective chemical potentials µ i (T ) of the stable hadrons in the partial chemical equilibrium case and for the temperature. Their knowledge is necessary for a determination of final state hadron abundances and spectra. We have restricted ourselves to the n B = 0 case in this work. In general, however, our approach allows for respecting the conservation of finite values for n B /s and n Q /s (while n S /s = 0) as relevant for heavy-ion collisions. Corresponding results for non-zero (although not too large) values of the associated chemical potentials will be reported in a forthcoming publication.