Constructing the Equation of State of QCD in a functional QCD based scheme

We construct the equation of state (EoS) of QCD based on the finite chemical potential information from the functional QCD approaches, with the assistance of the lattice QCD EoS. The obtained EoS is consistent with the up-to-date estimations of the QCD phase diagram, including a phase transition temperature at zero chemical potential of $T=155$ MeV, the curvature of the transition line $\kappa=0.016$ and also a critical end point at $(T,\mu_B)=(118, 600)$ MeV. In specific, the phase diagram mapping is achieved by incorporating the order parameters into the EoS, namely the dynamical quark mass for the chiral phase transition together with the Polyakov loop parameter for the deconfinement phase transition. We also implement the EoS in hydrodynamic simulations to compute the particle yields, ratios and collective flow, and find that our obtained EoS agrees well with the commonly used one based on the combination of lattice QCD simulation and hadron resonance gas model.


I. INTRODUCTION
One of the main goals of the relativistic heavy ion collision experiments is to search for the critical end point (CEP) of the QCD matter and to explore the thermodynamic properties of the strong interaction matter at finite temperature and chemical potential [1][2][3][4][5][6].To build the bridge between these measurements and the theoretical studies, the key element is the equation of state (EoS) of the QCD matter [7][8][9][10][11][12][13], since the EoS plays a crucial role as an input in hydrodynamic simulations for the evolution of the fireballs produced in heavy-ion collisions.
On the theoretical side, the QCD phase diagram has been widely studied via effective models [14][15][16][17][18][19][20][21][22][23][24], and the QCD approaches like lattice QCD simulation [25][26][27], functional QCD methods including Dyson-Schwinger equations [28][29][30][31][32][33][34] and functional Renormalization group approach [35][36][37].Among them, the functional QCD methods have delivered a comprehensive investigation on the phase diagram and especially, on the location of the CEP at large chemical potential.The current computations have given the phase transition line which is consistent with the lattice QCD simulation at small chemical potential and provided an estimation of the CEP at about µ CEP B ≈ 600 to 650 MeV [32][33][34][35].However, there is still no complete computation for the EoS which matches these up-to-date results for the phase transition line and is accessible for the hydrodynamics simulations.The estimated CEP is not even incorporated in the commonly applied EoS.
In this article, we then take advantage of different theoretical approaches to obtain an improved construction on the EoS.It mainly takes the functional QCD results about the phase transition line at finite density, which include the critical temperature, the curvature and the location of the CEP at large chemical potential, and also with the assistance of the lattice EoS at zero chemical potential.The constructed EoS is also required to satisfy the constraints from the current results of the thermodynamic quantities.All the elements in the construction are formalized analytically, making the further application of the EoS accessible and convenient.We also incorporate our EoS into hydrodynamic simulations to compute the experimental observables like particle yields, their ratios and collective flow.
The article is organized as follows: In Sec.II, we present the framework of constructing and parameterizing the EoS.In Sec.III, we present the results of EoS in the (µ, T ) plane, and also the Maxwell construction in the first order phase transition which stablizes the evolution.Then in Sec.IV, we present the results of particle yields and ratios together with the elliptic flow v 2 after incorporating with the hydrodynamics simulation.In Sec.V, we summarize the main results and make further discussions and outlook.

II. CONSTRUCTION OF THE EOS IN A FUNCTIONAL QCD-BASED SCHEME
Lattice QCD simulations have delivered solid computations at vanishing chemical potential [38][39][40].However, due to the sign problem, it is difficult for lattice QCD to reach real chemical potentials directly.The Taylor expansion approach [41][42][43][44] or alternative expansion scheme [45][46][47][48] are also limited in small chemical potential region, and the latest lattice results can only cover up to µ B /T ≤ 3.5.On the other hand, the large chemical potential region is accessible by continuum QCD methods, i.e. functional QCD approaches including Dyson-Schwinger equations (DSEs) [28,31,49,50] and the functional renormalization group (fRG) method [36,51].Therefore, one may combine the advantages of the two methods.In this work this is done by taking the following inputs: the lattice EoS data at zero chemical potential, the phase transition line data which is consistent between the lattice and also functional QCD results at small chemical potential, and also a predicted location for the critical end point from functional QCD calculations.
To calculate the EoS of the QCD matter at finite temperature T and quark chemical potentials µ = {µ q }, we apply the integral relation between the QCD pressure P and the quark number densities {n q } [52-54]: where P (T, 0) is from the lattice QCD input.Here it is determined from the parameterized formula of the lattice QCD's trace anomaly I(T ) = (ϵ − 3P )/T 4 at vanishing chemical potential in Ref. [39]: (2) with t = T /(200 MeV), h 0 = 0.1396, h 1 = −0.1800,h 2 = 0.0350, f 0 = 2.76, f 1 = 6.79, f 2 = −5.29,g 1 = −0.47,g 2 = 1.04 for 2+1 flavor.The pressure at vanishing chemical potential is then obtained with: Therefore, the key quantity for the EoS at finite chemical potential is the quark number density.It is related to the quark propagator S q , which is the central element in the functional QCD approach: with ω n the Matsubara frequencies and p the spacial momentum; the trace is taken over the color index and the Dirac structure.One can further take the vanishing momentum approximation for the inverse quark propagator as: with M q the dynamical quark mass and A 0 the gluon condensate [55,56] which relates to the Polyakov loop Φ: Here for simplicity we take the gluon condensate only with the Cartan generator τ 3 following the Ref. [56], with the eigenvalues ϕ fund = {±φ/2, 0}.From this, a simple analytical form is available for the number density, which has great advantage on the further applications.Such construction is the main purpose of this article.The quark number density n q in the (T, µ q ) plane can be expressed by the dynamical mass M q and additionally, the Polyakov loop parameter from Eq. ( 6), as follows [57]: x ± (k; T, µ q ) = exp [(E q (k; T, µ q ) ∓ µ q )/T ] , The next step is to quantitatively incorporate the phase diagram into the two order parameters M q (T, µ q ) and Φ(T, µ q ).First of all, it has been shown from the lattice and functional QCD studies, that the chiral phase transition line takes the parametrized form as: For the (2+1)-flavor case, the state-of-the-art results are T c (0) = 155 MeV, κ = 0.016 and |λ| ≲ 10 −4 [25-27, 33, 35].In addition, functional QCD approaches estimate that the CEP locates at about µ CEP B ≈ 600 to 650 MeV [32,33,35,58].This gives then a strong constraint on the order parameter of chiral phase evolution.
In the present work we simply take µ CEP B = 600 MeV.The chiral symmetry is related to the dynamical mass generation of quarks which have been continuously verified in functional QCD approaches and also effective QCD models.The quark mass at vanishing momentum can be then regarded as the order parameter of chiral phase transition [29,59].In the vacuum, the typical mass scale for the light-flavor quarks is M 0 = 350 MeV which is suggested by lattice QCD [60,61] and functional QCD [35,62] calculations.At finite temperature and chemical potential, one can take the Ising parametrization [63,64] in which a CEP is included: The (T, µ B ) dependence of the order parameter can be mapped to the Ising parameter (r, h), and the order parameter can be parameterized as: with the typical parameters β = 0.326, δ = 4.80, h 0 = 0.394 and h(θ) = θ(1 − 0.76201 θ 2 + 0.00804 θ 4 ) as given in Ref. [63], and M 0 is the normalization constant.The complete mapping procedure is (T, µ B ) → (R, θ) → (r, h).In order to match the phase transition line Eq.( 11) precisely, we propose a non-linear map as follows: with .0 • and α 2 = α 1 + 90 • in our case.The mapping function f PT is calibrated so that at h = 0, Eqs. ( 16) becomes exactly the parametric equations of the phase transition line in Eq. (11), which requires (17) In other words, here we consider h as the "distance" towards the phase transition line, and r being the projection coordinate on the phase transition line.The unconstrained parameters here are ω and ρ, which are currently set the same as those reported in Ref. [63].To sum up, the parameters for the dynamical mass are listed in Table I and the respective temperature and chemical potential dependence of dynamical quark mass M q is depicted in the upper panel of Fig. 1.With such a setup the temperature and chemical potential dependence of the chiral order parameter is consistent with the results from lattice QCD simulations and functional QCD methods.12) to (16).for the light-flavor quarks, including the position of the CEP.
While M q is the order parameter of chiral phase transition and is parametrized above, the Polyakov loop parameter Φ is usually connected to the deconfinement transition.Here we bypass the argument of the relation between Polyakov loop and deconfinement, and simply consider Φ as an order parameter for another phase transition in the gluon sector.
We apply the Polyakov loop data at zero µ q from Ref. [57], which is denoted as Φ(T, 0) = L(t) with t = T /T c (0).The following fit function is used here as: scheme proposed in [48,65,66], which suggests a temperature scaling behaviour of the thermodynamic functions: Here we choose T c (0) and κ for the deconfinement phase transition the same as those for the chiral phase transition as shown in Fig. 1.Note that the two phase transition lines are not necessarily identical which may induce some new phases like quarkyonic phase [67,68].Since the QCD phase transition is mainly characterized by the chiral phase transtion, deviation of the two phase transitions may only induce some delicate changes in the observables which require further investigations in the future.
In short, we set up a functional QCD framework for the EoS by mapping the up-to-date knowledge of the QCD phase structure to the order parameters.The free parameters here are those in the map functions, namely ω, ρ and M 0 .As we will demonstrate in the following, the obtained EoS is comparable with the lattice QCD pre-diction and the experimental observations quantitatively, which has a weak dependence on the free parameters.

III. NUMERICAL RESULTS OF THE EOS WITH FIRST ORDER PHASE TRANSITION
Next, we consider the (2 + 1)-flavor case, the baryon, electric charge and strangeness chemical potentials (µ B , µ Q , µ S ) are associated with the quark chemical potentials as: We first consider the 3 flavor degenerate case with 3 which is the conventional case often applied in the hydrodynamics simulations.With the obtained pressure P (T, µ B ) and number density n q (T, µ B ) and using thermodynamic relations, one can compute the entropy density s = ∂P/∂T , the energy density ε = T s − P + µ B n B and so on.The speed of sound squared c 2 s is then given by c 2 s = ∂P ∂ε .We show thus the 3-dimensional (3D) plot in terms of temperature and chemical potential for the pressure, energy density, number density and speed of sound in Fig. 2. The current results are consistent with the phase transition line obtained by lattice QCD simulations and functional QCD methods with a CEP at (T, µ B ) = (118, 600) MeV.Now one may take a closer look at the EoS at different chemical potentials as depicted in Fig. 3.For small chemical potentials, the speed of sound is a smooth function with a minimum around c 2 s ∼ 0.12 at phase transition point.As the chemical potential becomes larger, the speed of sound becomes more oscillated near the transition temperature.At the CEP and the first order phase FIG. 2. Calculated 3D plots for the pressure, number density, entropy density and energy density in terms of the temperature and chemical potential, normalized to a dimensionless form by the temperature T with the respective power.transition region, the speed of sound at phase transition point becomes zero.For the first order phase transition, since one has dP = 0 and a finite dε during the phase transition, the curve becomes discontinuous for the two phases.The speed of sound in the vicinity of phase transition point is still finite, and only a few points reach to zero drastically.More interestingly, as the chemical potential increases, there appears a peak nearly above the phase transition temperature.The maximum value of the peak gradually grows and saturates to the conformal limit c 2 s = 1/3, which implies a new feature of the QCD EoS at high densities that may shed light on the studies of neutron stars [69][70][71].
In the ideal case of the first order phase transition, there exists discontinuity in the number density, entropy density and also energy density.The discontinuity would be affected by the non-equilibrium effects in the dynamical evolution and will cause problems in the hydrodynamics simulation.Here, we consider the Maxwell construction to fill the discontinuity with simply a linear transition of the thermodynamic functions from one phase to the other.This construction leads to a vanishing speed of sound during the first order phase transition.To consider the spinodal decomposition, one may need a construction described in Ref. [63], and the speed of sound will be negative in the instable region.In Fig. 4, we show the construction for the pressure as the function of energy density.Using this construction, one can convert the thermodynamic quantities from the (T, µ B ) dependence into the (ϵ, n B ) dependence, which is then accessible for the hydrodynamic simulations.It also needs to mention that based on the current feature of EoS, a large value of c 2 s near conformal limit does not necessarily mean that there is no first order phase transition in neutron star as indicated in Ref. [72], on the contrary, it might be an important feature when the QCD matter approaches the first order phase transition region at large chemical potential.The EoS data as a function of T and µ B together with the MUSIC input format for both B and BS cases are available in the github (fQCD-EoS-PhaseDiagramMap).
At last, we check the conserved charge conditions satisfied in the heavy-ion collisions.Consider the 3 light flavors f = u, d, s, the conditions are: with conserved charges n B,Q,S which stand for the baryon, electric and strangeness densities, respectively, and r the charge-to-mass ratio of the collided nuclei, e.g.r Au+Au ≈ 0.4.
In this sense the 3-flavor-degenerate scenario is equivalent to µ Q = µ S = 0 in Eqs. ( 21) to (23), i.e.only the baryon chemical potential is considered.On the other hand, for the (2+1)-flavor scenario Eq. ( 24) is equivalent to: Within our framework Eq. ( 25) requires µ s = 0. Therefore, for the charge conserved case the remaining task is to consistently solve Eqs. ( 21), ( 22), ( 25) and (26).The parameters in Table I is adapted for both u and d quark densities in Eq. ( 26).we show the calculated isentropic trajectories with s/n B = const.in Fig. 5, for the 3-flavor-degenerate case which is denoted as B (baryon), and the charge-conserved case with the constraint of the conserved charge conditions which is denoted as BS (baryon-strange).The obtained isentropic trajectories are in good agreement of the previous studies.It shows that the strangeness neutrality pushes the trajectories into higher chemical potential region, which may have a impact on the initial conditions of heavy ion collisions [73][74][75][76].

IV. APPLICATION IN HYDRODYNAMIC SIMULATION
To evaluate the applicability of the constructed fQCD EoS in realistic dynamical evolution, we carry out the hydrodynamic simulation with the EoS as input.In our hydrodynamic simulations, we implement the hydrodynamic model MUSIC [77,78] with the initial condition generated by AMPT [79], where the initial energy/baryon number distribution are extracted at constant proper time τ 0 .To compare the difference between various EoSs, we calculate the observables of thermal particles without hadronic scattering or decay effects.See Ref. [80] for more details.Here we do the calculation with ideal hydro and neglect the initial flow, which do not affect the final particle yields obviously.Fig. 6 shows the identified particle yields and ratios calculated by AMPT + MUSIC in 0-5% Au+Au collisions at √ s N N = 39 GeV.As a comparison, we also show the results calculated by a crossover equation of state NEOS-B, which is based on lattice QCD simulation and hadron resonance gas model [81].Note in most central Au+Au collisions at 39 GeV, the corresponding chemical freeze-out (T, µ B ) is around (156, 160) MeV, where fQCD-B and NEOS show no obvious difference.When considering the charge conservation in fQCD-BS, the anti-baryon/baryon ratio decreases due to larger baryon chemical potential, which also corresponds to the isentropic trajectories in Fig. 5.The ratio for hyperon like Λ and Ξ also show similar behavior.Note that here the fQCD EoS-BS changes the anti-baryon/baryon ration, which may be further improved by considering the flavor mixing effect in the future.In addition, Fig. 7 shows the corresponding elliptic flow v 2 of π + and proton calculated by pure hydrodynamics.It shows evidently that, at the crossover region, T h e r m a l p a r t i c l e d N / d y ( y = 0 ) T h e r m a l p a r t i c l e r a t i o A u A u , 3 9 G e V , 0 -5 % FIG. 6.
Identified thermal particle yields (upper panel) and their ratios (lower panel) in 0-5% Au+Au collisions at √ sNN = 39 GeV.The scatters represent the model calculations with NEoS-B and fQCD equation of states on freeze-out surface, respectively.changing EoS will not make obvious difference to the bulk evolution, which is also consistent with the result of Fig. 5.
To study the stability of the fQCD EoS near the first order phase transition, we run the hydrodynamics at higher baryon chemical potential region without freezeout.Fig. 8 shows the obtained averaged hydrodynamic trajectories at such baryon-rich region, where we take the AMPT initial profile at √ s N N = 7.7 GeV but re-scaled the initial net baryon density distribution by factors to match s/n B ≈ 13 and s/n B ≈ 8 in initial conditions, respectively.In both cases, hydrodynamics start with same initial energy profile and can run stably for more than 10 fm/c.Note for the case s/n B ≈ 8, the theo- retical isentropic line will experience a first order phase transition with decreasing temperature.While in hydrodynamic simulation, the speed of sound will turn to zero at phase transition, so it cannot evolve across the phase transition line but go along with it instead.How the system evolves after first phase transition is out of the reach of current hydrodynamic model, a more realistic dynamical description applying the constructed EoS with CEP and first order phase transition will be investigated in the future.

V. SUMMARY
We proposed an improved construction on the QCD EoS in a functional QCD based scheme, which takes input from the lattice EoS at zero density, and moreover utilizes the knowledge of the QCD phase diagram at finite density from functional QCD approaches, i.e. the phase transition line and the location of the CEP.By implementing the zero-momentum approximation, the quark number density is expressed analytically, which ensures then that the EoS can be analytically calculated and is convenient for further applications.
We computed the thermodynamic quantities such as the pressure, the energy density and the entropy density, and eventually the speed of sound as a function of temperature and chemical potential.For small chemical potentials, the obtained minimum around c 2 s ∼ 0.12 at the phase transition point is consistent with the results from lattice QCD simulation.As the chemical potential increases, the speed of sound drops drastically to zero at the phase transition point.Moreover, the speed of sound at large chemical potential shows a conformal limit behavior right after the first order phase transition, which may have some impacts on neutron star properties.
We applied our constructed EoS to get the isentropic trajectories.Our results are in a good agreement with those from lattice QCD simulations.Besides, the current estimation gives that the evolution goes into first order phase transition region at s/n B ∼ 10.This provides a benchmark for the QCD chiral phase transition in heavy ion collisions and the possible new phenomena such as the non-monotonicity of the triton yield ratio.The further investigation either resolve the conflicts between theory and experiment, or reveal some new features of QCD.We also put the EoS into the hydrodynamic simulations.For Au+Au collisions at √ s N N = 39 GeV, our results on particle yields, their ratios and also the elliptic flow of pion and proton are in good agreement with the results from a commonly used EoS.Moreover, we showed that with the current EoS, it is possible for the hydrodynamic simulations to reach the first order phase transition region.A more realistic dynamical description of the first order phase transition will be investigated in the future.

with l 1 1 ΦFIG. 1 .
FIG. 1. Density plots for the two order parameters Mq and Φ as functions of temperature T and chemical potential µq for light-flavor quarks.

FIG. 3 .FIG. 4 .
FIG. 3. Calculated speed of sound squared as a function of temperature T at several values of the baryon chemical potential µ B .

FIG. 8 .
FIG.8.Hydrodynamic trajectories at high baryon number density region.The scatters shows the hydro evolution with different initial s/n, where they are averaged in each time step with energy density weight.The blue "star" shows the critical end point locate at (T, µB) = (118, 600) MeV.

TABLE I .
The mapping parameters of the dynamical mass in Eqs. (