Isentropic evolution of the matter in heavy-ion collisions and the search for the critical endpoint

We study the isentropic evolution of the matter produced in relativistic heavy-ion collisions for various values of the entropy-per-baryon ratio of interest for the ongoing and future experimental searches for the critical endpoint (CEP) in the QCD phase diagram: these includes the current Beam-Energy-Scan (BES) program at RHIC and the fixed-target collisions foreseen for the near future at various facilities. We describe the hot-dense matter through two different effective Lagrangians: the PNJL (Polyakov-Nambu-Jona-Lasinio) and the PQM (Polyakov-quark-meson) models. We focus on quantities expected to have a direct experimental relevance: the speed of sound, responsible for the collective acceleration of the fireball, and the generalized susceptibilities, connected to the cumulants of the distributions of conserved charges. In principle, they should affect the momentum spectra and the event-by-event fluctuations of the yields of identified particles. Taking realistic values for the initial temperature and the entropy-per-baryon ratio we study the temporal evolution of the above quantities looking for differences along isentropic trajectories covering different regions of the QCD phase diagram, passing far or close to the CEP or even intersecting the first-order critical line.


Introduction
The goal of relativistic heavy-ion collisions is the exploration of the QCD phase diagram, looking for signatures of the transition from a system of colour-neutral hadrons to a deconfined plasma of quarks and gluons (QGP). Actually, in current ongoing experiments at RHIC and at the LHC, the transition occurs in the opposite direction: at very high energy, in the collisions of two nuclei, the amount of stopped baryonic matter is negligible, quarks and gluons are "newly" produced particles arising from the strong colour fields in the overlapping region. Quarks and gluons form a thermalized plasma undergoing an almost adiabatic expansion during which the latter cools down until reaching a temperature at which colour-singlet hadrons become again the active degrees of freedom.
First-principle lattice-QCD simulations show that, at vanishing baryon density (i.e. at baryo-chemical potential µ B = 0), the transition connecting the partonic and hadronic phases is actually a smooth crossover [1]. This is the regime of relevance for the nuclear collisions at the LHC ( √ s NN = 2.76 and 5.02 TeV) and at the highest center-of-mass energy at RHIC ( √ s NN = 200 GeV) and this corresponds also to the regime at which the QCD transition occurred during the thermal history of the uni-verse, around 1 µs after the Big Bang, when the temperature reached a value around 150-160 MeV [2]. Unfortunately, due to the sign problem which prevents a Monte-Carlo sampling of the gauge-field configurations, lattice-QCD simulations cannot provide definite answers on the QCD thermodynamics and phase structure at finite baryon density, except for sufficiently small values of µ B /T where, for instance, one can perform a Taylor expansion around µ B = 0 [3,4]. Various effective models were then employed (see e.g. [5,6,7]), suggesting that at large µ B the transition should become of first order. Before turning into a crossover, such a first-order line in the (µ B , T )-plane should end with a critical endpoint (CEP) in which the transition is of second order, characterized by an infinite correlation length.
The study of the QCD phase-diagram in the region of non-vanishing baryon density (for a recent review see for instance [8]) presents then several challenges. First of all one should provide a non-ambiguous, quantitative theoretical answer to the question concerning the existence and location of a CEP; so far lattice-QCD simulations tend to disfavour its presence for µ B /T < ∼ 2 [4]. Secondly, one should suggest specific signatures of the QCD critical point to look for in the experimental data [9,10]. Finally, one should perform experiments able to explore the QCD matter in such a regime of high baryon-density. This can be achieved through the collisions of heavy nuclei at lower center-of-mass energies, characterized by a larger stopping of the incoming nucleons. For this purpose, experiments were performed in the past at AGS [11] and SPS [12], the Beam-Energy-Scan (BES) is currently ongoing at RHIC [13,14], a rich physics program is scheduled for the present and the near future at SPS [15] and new infrastructures like NICA [16] and FAIR [17] are under construction. One of the difficulties in identifying unambiguous experimental signatures of the presence of a CEP and/or of the occurrence of a first-order transition in the hot QCD matter produced in relativistic heavy-ion collisions is that -at variance with condensed-matter experiments -one does not have the possibility to explore the phase diagram in a systematic way, tuning appropriate control parameters and performing measurements during the whole evolution of the system (except for the case of rare penetrating probes like photons and dileptons). One can only change the center-of-mass energy of the collisions and the nuclear species and possibly select events of different centrality; furthermore, one can only measure hadrons after their chemical and kinetic decoupling from the expanding fireball.
In heavy-ion collisions one produces a system whichneglecting dissipative effects due to viscosity, heat conduction and charge diffusion -undergoes an approximate isentropic expansion moved by pressure gradients along trajectories of constant entropy per baryon (s/n B = const), the higher the center-of-mass energy of the collision the higher the s/n B ratio. Hence, in performing theoretical calculations of experimental relevance, the quantities of physical interest must be evaluated along the above trajectories. This is what we plan to do in our paper, focusing in particular on two kind of thermodynamic quantities: the speed of sound and the generalized quark-number susceptibilities. Both of them have a deep experimental relevance. The squared speed of sound c 2 s governs the response of the system to the initial energy-density gradients, leading to the collective acceleration of the fireball. If, during its evolution, the system undergoes a first-order transition during the mixed-phase the value of the speed of sound should be very low and this should affect the transverse-momentum distributions of the hadrons produced in the collision. The generalized susceptibilities, on the other hand, are related to the fluctuations of conserved charges (baryon number, electric charge and strangeness) which can be experimentally accessed through event-byevent measurements of the unbalance of protons and antiprotons [18], opposite-charge particles [19] and, more recently, K + and K − mesons [20]. They are expected to display huge oscillations and rapid changes of sign in the vicinity of a CEP, and one hopes that this will leave observables effects in the event-by-event measurements of identified particles produced at the chemical freeze-out, if the latter occurs sufficiently close to the CEP.
In our paper we evaluate the above quantities starting from two effective Lagrangians which display the same chiral-symmetry breaking/restoration pattern of QCD and include a simple modelling of quark confinement: the PNJL (Polyakov-Nambu-Jona-Lasinio) model [21] and the PQM (Polyakov-Quark-Meson) model [22]. We consider isentropic trajectories of relevance for the BES undergoing at RHIC and for future experiments at SPS, NICA and FAIR. The corresponding values of s/n B are estimated starting from the data provided by the STAR collaboration in [14,23]. We study the values taken by the speed of sound and by the generalized susceptibilities during the isentropic evolution of the system until the estimated freeze-out point, looking for differences among cases in which the QCD transition occurs via a smooth crossover, a first-order transition or crossing a possible CEP.
Our paper is organized as follows. In Sec. 2 we present the effective models employed for our calculations. In Sec. 3 we display our results for the phase diagram of the two models and for the evolution of the speed of sound and of the generalized quark susceptibilities along very different isentropic trajectories. Finally, in Sec. 4 we discuss the obtained results and their possible experimental relevance, providing also some perspective for future developments.

The PNJL and PQM models
In this section we present the two effective Lagrangians used in our study: the PNJL and the PQM model. The two models share some similarities, in particular in the effective implementation of quark confinement, but they differ in the way in which the interaction among quarks is described: via a Fermi-like four-fermion vertex in the PNJL model, via the exchange of scalar/pseudoscalar mesons representing bosonized quark fields in the PQM model. Both models have been already described at length in the literature and for more details we refer the interested reader to the original publications on the subject [21,22].

The Polyakov field and the effective implementation of confinement
In QCD -for infinite quark mass -the order parameter of the confinement/deconfinement phase transition is the Polyakov loop, expressed in terms of the Euclidean temporal component of the gauge (1) Its expectation value is related via the equation to the change in free-energy occurring adding an isolated heavy quark into the system: it vanishes in the confined phase, where adding an isolated colour charge requires an infinite amount of energy; it is non-zero in the deconfined phase, where this suppression of coloured states is absent. From the mathematical point of view Φ = 0 leads to the spontaneous breaking of the global Z Nc symmetry, related to the transformations which multiplies all the Polyakov lines L x by the same phase factor z k . In the pure SU (N c ) gauge theory this is an exact symmetry of the Yang Mills action. Light dynamical quarks introduce an explicit breaking of the Z Nc symmetry (N c = 3 in QCD), nevertheless the Polyakov loop, in light of its physical meaning, remains a useful quantity to identify the different phases predicted by the theory. Both in the PNJL and in the PQM model the Polyakov loop Φ (andΦ for the antiquark sector) plays the role of a constant background field whose value is determined looking for the stationary points of the thermodynamic potential and whose role is to suppress the contribution of free quarks to the thermodynamic quantities. In order to ensure the occurrence of a phase transition, an effective temperature-dependent potential for the Φ andΦ fields is introduced into the Lagrangian, with the constraint of respecting the Z 3 symmetry of the theory. In this paper we choose it of the following polynomial form: and the coefficients of the potential are fitted in order to reproduce the lattice data for pure Yang-Mills Theory [21]: their values are given Table 1

The PNJL and PQM Lagrangians
The PNJL and PQM Lagrangians can be written as L = L 0 + ∆L, where L 0 is the piece common to both models (including the quark kinetic term, the coupling with the background Polyakov field and the Polyakov potential), while ∆L describes the quark interaction, either through a contact term (in the PNJL model) or mediated by light scalar/pseudoscalar mesons (PQM model). We start introducing the common piece of the Lagrangian, which reads In the Dirac term the covariant derivative is defined as is the SU c (3) gauge field and λ a , a = 1, .., 8 are the Gell-Mann matrices of the SU c (3) group. The symbol µ f indicates the chemical-potential matrix associated to the various quark flavours, µ f ≡ diag(µ u , µ d , µ s ), which are conserved by strong interactions. Notice that the structure of the Lagrangian entails that the background gauge field enters into the quark propagator as a shift of the chemical potential, µ −→ µ + A 0 ≡ µ + iA E 4 . Finally, the term U Pol is the Polyakovloop potential given in Eq. (4).
It is useful to remind the relation between the quark chemical-potentials and the ones associated to the conserved QCD charges at the hadronic level, which are the ones of actual experimental relevance: Here µ B , µ Q , µ S are the baryon, electric-charge and strangeness chemical potentials. We now discuss the interaction terms of the Lagrangian, in which the two models differ.
3.67/Λ 2 −12.36/Λ 5 602.3 5.5 140.7 Table 2. Parameters of the NJL model We start from the PNJL model, for which we have: In the above m 0 is the diagonal mass matrix for the quarks, m 0 ≡ diag(m 0u , m 0d , m 0s ), λ a (a = 1, 2, ..., 8) are the 3 × 3 Gell-Mann matrices of the SU f (3) group and λ 0 = 2/31. G and K are the coupling constant for the 4fermion vertex and 6-fermion vertex, with dimensions [E] −2 and [E] −5 , respectively. The six-fermion term is know as 't-Hooft determinant and is introduced to explicitly break the axial U A (1) symmetry, which is not a symmetry of QCD at the quantum level, as confirmed by the high mass of the pseudoscalar η meson. Due to the above interaction terms the PNJL model in 4 space-time dimensions is not renormalizable and it has to be considered just a low-energy effective theory. To regularize the divergent integrals one has to introduce a momentum cutoff Λ, representing a further free parameter of the model together with the couplings G and K and the quark masses: they are all summarized in Table 2. We work in the isospinsymmetric case, with m 0u = m 0d ≡ m 0 , setting also Considering now the PQM model, the terms of the Lagrangian describing the scalar and pseudoscalar mesons and their interaction with the quark field read where φ combines the scalar and pseudoscalar meson fields g denotes the Yukawa coupling between quarks and mesons, k is the coefficient of the U A (1) symmetry breaking term and H = T a h a . The coefficients k, h a , m 2 , λ 1 and λ 2 are the parameters of the mesonic Mexican-hat potential [24,25,26] which describes the spontaneous breaking of chiral symmetry, giving the meson fields their expectation valuesσ a . The explicit breaking of chiral symmetry in the PQM model arises from the tilt of the mesonic potential provided by the term H(φ + φ † ), playing the role of the −ψ m 0 ψ term in the PNJL Lagrangian. The quantities used to fix the parameters of the PQM model are summarized in Table 3.

The mean-field approximation and the thermodynamic potential
In this section we display the equations allowing one to study the thermodynamics of the PNJL and PQM models in the mean field approximation, at the basis of the numerical results presented in Sec. 3. Within this setup a system of strongly-coupled quarks is described as a collection of non-interacting quasi-particles, endowed with effective masses obtained through the minimization of the resulting thermodynamic potential.
In the two models the effective quark mass has a different origin. In the PQM case it arises from the Yukawa coupling with the meson fields, which develop a non-vanishing expectation value, leading to m = gσ and m s = √ 2gσ s for the light and strange quarks, respectively [27]. In the PNJL case it arises from the quark self-interaction and is obtained linearizing the Lagrangian given in Eq. (8) with respect to the fluctuations of the composite operator and hence the effective quark masses (we assume exact isospin symmetry, setting m u = m d ≡ m ) turn out to be expressed in terms of the chiral condensates ϕ f ≡ ψ f ψ f of the various quark flavours Getting the quark masses requires the self-consistent solution of the above system of coupled gap equations, where the mass of quarks of flavour i depends on all chiral condensates Here E i p = p 2 + m 2 i is the energy of the dressed quark of flavour i, and the modified Fermi distributions are given by and Also the thermodynamic potential can be conveniently written as Ω ≡ Ω 0 + ∆Ω, where Ω 0 is the part common to both models, while ∆Ω includes the model-dependent terms. At the mean field level one has having defined and (18) Concerning the model-dependent contribution to the thermodynamic potential, in the PNJL case one has: For the PQM model the contribution from the expectation value of the meson fields reads: At the mean field level the thermodynamic potential is a function of the effective quark masses and of the expectation value of the Polyakov fields, which have to be self-consistently determined requiring Ω MF to be stationary under variation of the above quantities, i.e.
The first two conditions lead to the mass-gap equations (13), which in the PNJL model were independently obtained expressing the chiral condensates ϕ i in terms of the Hartree quark propagators. Equivalently, the dependence on the effective quark masses m and m s can be traded for the one on the chiral condensates ϕ and ϕ s or on the expectation values of the meson fields σ and σ s . Notice that in the present MF approximation, replacing the meson fields by their expectation values, the associated kinetic contribution in the PQM Lagrangian in Eq. (9) vanishes. A final comment on the UV behaviour of the two models is in order. The second term in the RHS of Eq. (16) contains the sum of the (negative) zero-point energies of the various fermionic modes. It is divergent and is regularized by the UV cutoff Λ. In the PNJL case the results depend explicitly on this cutoff, which -as already mentionedhas to be viewed as a parameter of the model, fixed by matching some zero-temperature/density observables. On the other hand, the PQM model is renormalizable and in principle one could cancel the dependence on Λ of the vacuum fluctuations with the dependence on the cutoff of the parameters of the meson potential in Eq. (20), so that the final physical results do not depend on the choice of the UV regulator. Here however, in order to perform a consistent comparison, we simply regularize the UV divegence through the same ultraviolet cutoff Λ employed in the PNJL model.

Results
Our results focus on the phase-diagram of the PNJL and PQM models, on the nature of the chiral-deconfinement transition, on the location of the critical endpoint (CEP), on the speed of sound and on the generalized susceptibilities, the latter being of interest since associated to the thermal fluctuations of conserved charges. All our numerical calculations are performed along isentropic trajectories. The comparison of the angular distributions of soft identified hadrons with the results of hydrodynamic calculations suggests in fact a very low value of the viscosity to entropy-density ratio, compatible with the conjectured universal lower bound η/s = 1/(4π) predicted by the gauge-gravity duality. Dissipative effects are then small in the fireball produced in heavy-ion collisions, which evolves along trajectories of constant entropy per baryon: both the entropy and the baryon density get dilute due to the expansion of the system (which has a positive pressure with respect to the surrounding vacuum), but their ratio s/n B remains constant within each fluid-cell.
In order to provide results of phenomenological relevance we need to estimate both the entropy-per-baryon ratio S/B and the initial entropy density of the system arising from the heavy-ion collisions at the various nucleonnucleon center-of-mass energies explored in the BES at RHIC, from 7.7 GeV to 200 GeV. We start estimating the initial entropy density s 0 , assuming that its value is proportional to the measured rapidity density of charged particles dN ch /dη. For Au-Au collisions at √ s NN = 200 GeV the value s 0 = 84 fm −3 , once inserted in the initial condition of hydrodynamic calculations, was shown to satisfactory reproduce soft-hadron distributions [28] and, later on, it was widely employed in the literature, e.g. [29,30]. Taking τ 0 = 1 fm/c as an estimate of the initial thermalization time and integrating over the transverse plane the profile provided by a Glauber calculation one gets dS/dy ≈ 4700 for the entropy per unit rapidity in central Au-Au collisions; taking into account that, for a pion gas around the chemical freeze-out temperature, one has S > ∼ 4N , this compares well with the observed rapidity density of charged particles [31]. The initial entropy density s 0 at lower center-of-mass energies is obtained rescaling the estimate at √ s NN = 200 GeV according to the lower values of dN ch /dy [14,23]. Also the S/B ratio at the various center-of-mass energies is estimated from the yields of identified hadrons -π ± , K ± , p/p -quoted in Refs. [14,23], still assuming that each particle carries about 4 unit of entropy. Our results for s 0 and S/B are collected in Table 4, where we also quote the values of the kinetic freezeout temperatures obtained in [14,23] through a blast-wave fit of the transverse-momentum distributions of identified hadrons.

Phase diagram and isentropic trajectories
We start our analysis with the study of the phase diagram of the two models and of the evolution of the system along isentropic trajectories corresponding to values of s/n B of  Table 4. Each curve starts at a temperature T0 corresponding to the values of s0 given in Table 4 and is plotted down to the corresponding kinetic freeze-out temperature T fo . We also plot two isentropic trajectories so far not covered by the BES and entering into the first-order region. Fig. 2. The same as in Fig. 1 Table 4. Estimate of the initial entropy density and of the entropy per baryon in Au-Au collisions at different center-ofmass energies. We also quote the kinetic freeze-out temperature obtained in Refs. [14,23] through a blast-wave-fit.
interest for the ongoing BES at RHIC. We display our results both in the µ q −T and in the n q −T planes, the last case allowing one to get a deeper physical insight when the system, during its evolution, meets a first-order phase transition. Our findings are shown in Figs. 1 and 2. There is a wide region in the phase diagram in which the transition from the chirally-symmetric, deconfined phase to the chirally-broken, confined one is actually a crossover. In this case there is no well-defined location of the transition associated to an unambiguous order parameter, but there are several quantities displaying peaks in a quite narrow range of temperatures and chemical potentials. We decide to identify the crossover between the two phases -shown as a black, dotted line in Figs. (1) and (2) -through the inflection point of the effective mass of the light quarks. We also display the isentropic trajectories followed by the system corresponding to the s/n B values quoted in Table 4; for each case, they are plotted starting from s 0 down to the estimated kinetic freeze-out temperature T fo kin . In both the PNJL and the PQM models, for the currently accessible values of s/n B , the transition of the system between the two phases occurs in the crossover region.
The crossover line ends at a critical endpoint, beyond which the transition becomes of first order, displayed by black solid lines in the figures. The location (µ CEP B , T CEP ) of the CEP is found to be (875,121) MeV and (903,118) MeV in the PNJL and PQM model. This can be compared with the results found in independent implementations of the PNJL model [32]. The corresponding critical isentropes passing through the CEP is given by values of the entropy per baryon s/n B = 7.02 and s/n B = 6.16, respectively. We decide then to show also a few isentropes in which the system, during its evolution, either passes very close to the CEP or meets the first-order line. This region is shown in finer detail in Fig. 3. In this last case, along the critical line, for a given T and µ q the thermodynamic potential has two degenerate minima. The single critical line in the µ q − T plane actually becomes the boundary of an extended region of phase coexistence in the n q −T plane. In the coexistence region no stable homogeneous phase can exist, but a fraction α of the volume is occupied by the chirally-restored phase and a fraction (1 − α) is occupied by the chirally-broken one. They are characterized by the same pressure, temperature and chemical potential, expressing the mechanical, thermal and chemical equilibrium between the two phases. For each T and µ q the value of α, which is determined through a Maxwell construction, depends on the history and kind of evolution of the system. Usually in thermodynamics one considers phase transitions occurring along isothermal lines, however in heavy-ion collisions there is no thermal bath with which the fireball is in contact. The system follows then an isentropic expansion at fixed s/n q (n q = n B /3) and one has where s Q/H and n Q/H are the entropy and quark-number density in the chirally restored/broken phase. One has α = 1 when the isentrope crosses the high-density branch of the critical line and α = 0 when it crosses the low-density branch. Hence, a given value of T and µ q does not fix α: one has to specify also the isentrope followed by the system. Notice that, considering this kind of evolution, one is implicitly assuming that the nucleation rate of bubbles of the low-density, chirally-broken phase is larger than the expansion rate of the fireball and we are not addressing the case of super-heating/cooling, which occurs if the system remains in a metastable minimum. Establishing whether this is a realistic assumption of whether the transition occurs via a different dynamics (spinodal decomposition) would need a deeper study of the rate of bubble nucleation, requiring in particular the evaluation of the surface tension of the interface between the two phases [33,34,35,36], which is out of the scope of the present paper. In any case, in Figs. 1 and 2, we display the regions in which in the two models a homogeneous metastable phase may exist, extending from the continuous first-order line to the dashed spinodal lines. The phase diagram plotted in terms of the quark density rather then of the chemical potential reveals how, in both models, the first-order coexistence region extends down to the origin. In the low-temperature regime the behaviour of both models is then unphysical, since they do not leave room for the existence of a self-bound homogeneous nuclear-matter phase, which we know to exist, but whose experimental density n B ≈ 0.16 fm −3 would lie here in the coexistence region. Furthermore there is no room for the liquid-gas phase transition of nuclear matter, which is a characteristic feature of strong interactions in the low-temperature (T < ∼ 20 MeV), high-density regime for low values of the entropy per baryon [37,38,39]. This must be viewed as a shortcoming of the models due to the pure scalar/pseudoscalar interaction and to the mean field approximation. We expect that the inclusion of a vector interaction and of hadrons as dynamical degrees of freedom in the confined, chirally-broken phase should improve the description of this region of the phase diagram.
In order to assess the role of the Polyakov field in fixing the location of the CEP in the phase diagram, in Fig. 4 we compare the PNJL vs NJL (left panel) and PQM vs QM models (right panel). In both cases, coupling the quark with the Polyakov field, the CEP turns out to move to a higher temperature and a lower baryo-chemical potential.
Finally, in Fig. 5, we show a few isothermal curves for the two models, plotting the pressure as a function of the specific volume 1/n B . We choose values of T either above or below the critical value T CEP . In this last case, in which a first-order transition occurs, we show also the Maxwell construction. Notice that the isothermal curves with T < T CEP display two stationary points. Between the two stationary point the pressure is a decreasing function of the density and the system is thus unstable. This part of the curves actually corresponds to effective quark masses arising from the maximum of the thermodynamic potential; this is also the region in which a single homogeneous phase can not exist.

The speed of sound
When the high temperature of a system allows the continuous creation/annihilation of particle-antiparticle pairs with m T the mass density in no longer a meaningful concept and the Newtonian definition of speed of sound has to be accordingly generalized. One can show that the relativistic squared speed-of-sound c 2 s is then given by the following derivative of the pressure with respect to the energy density at constant entropy per particle (baryon in the case of strong interactions): The speed of sound maps then a density gradient into a pressure gradient, which -if the evolution of the system can be described by hydrodynamics -is responsible  for the acceleration of the fluid. Having already evaluated the various isentropic trajectories at s/n B = const, the above quantity is simply obtained via a numerical finite-difference method. Results for c 2 s corresponding to the values of s/n B considered in this work are shown in Figs. 6 and 7, plotted as a function of temperature T and time τ . In mapping T into τ we assume a simple Bjorken-like inviscid longitudinal expansion, for which one has sτ = s 0 τ 0 . The curves are plotted down to the kinetic freeze-out temperatures T fo kin provided by the experimental blast-wave fits and quoted in Table 4. Since -in the meanfield approximation -as active degrees of freedom in the low-temperature/density phase both models have dressed quarks suppressed by their large effective mass and by their coupling with the Polyakov field, a kinetic freeze-out temperature T fo kin ≈ 120 MeV corresponds to a very small value of the entropy density, much smaller than the one of a hadron-resonance gas; hence, in order to reach the experimental T fo kin , we have to follow the evolution of the system up to unrealistically large values of time τ . Having a more realistic description of the chirally-broken phase within the two models would require to include as dynamical degrees of freedom the set of scalar/pseudoscalar mesons arising from the quark-antiquark interaction contained in the Lagrangian. We leave this issue for future work.
In Figs. 6 and 7 one can appreciate the very different behaviour of the speed of sound depending whether, during the isentropic evolution, the transition from the chirally restored to the broken phase occurs via a smooth crossover (curves with s/n B ≥ 17.5) or via a first-order change of state (curves with s/n B = 2, 5). For the s/n B values of experimental relevance in both models the transition occurs in the crossover region and the speed of  Fig. 1 plotted as a function of temperature (left panel) and time τ (right panel). The different behaviour between a first-order transition (s/nB = 2, 5) and a smooth crossover is clearly visible. Fig. 7. The same as in Fig. 6, but for the PQM model. sound displays a rapid but continuous drop for temperatures around the inflection point of the light chiral condensate. Notice that there is no qualitative difference in the behaviour of the curves corresponding to the various center-of-mass energies explored at RHIC, although -at a quantitative level -for all values of time the speed of sound for low s/n B is significantly lower than the one for large s/n B : this should have an impact on the final momentum distributions of soft hadrons.
On the other hand if, during its evolution, the system meets a first-order transition the speed of sound suddenly drops (almost) to zero. Such a soft equation of state (EoS) should translate into a very small acceleration of the fluid during this stage, leaving its signatures in the final momentum distributions of the hadrons decoupling from the fireball.
More quantitative considerations would require inserting the above EoS into a full hydrodynamic code including also the evolution of baryon density [40], allowing one to simulate the expansion of the fireball and the continuous decoupling of light hadrons.

Generalized susceptibilities
In this section we evaluate higher-order cumulants of the distributions of conserved charges (in particular, baryon number) provided by the PNJL and the PQM models. For the sake of consistency, here we remind the basic definitions. In statistical mechanics cumulants are expressed in terms of derivatives of the pressure with respect to the chemical potential, i.e [41].
whereN is the conserved charge and µ the associated chemical potential. In heavy-ion collisions the relevant conserved charges are baryon number, strangeness and electric charge, however here we simply focus on baryon-number fluctuations. The first few cumulants of a distribution are given by (δN ≡N − N ) (25) and allow one to define its mean M , its variance σ 2 , its skewness S and kurtosis κ (26) In particular, skewness and kurtosis of a distribution quantify how asymmetric and peaked/broad the latter is with respect to a Gaussian, for which N n c = 0 for n ≥ 3 and hence S = κ = 0.
Why is the study of higher order cumulants of interest for heavy-ion collisions? A first motivation comes from the search of the CEP in the QCD phase-diagram, where -if present -the transition is of second order and hence characterized by an infinite correlation length of the order parameter, in this case the chiral condensate (ψψ). One finds that higher order cumulants display a stronger sensitivity to such a quantity [9] and this is of relevance for experimental studies, in which the finite size and the expansion rate of the produced medium -far from the thermodynamic limit in which a phase transition is rigorously defined -would prevent one from observing any actual divergence of a correlation length. Secondly, the experimental measurement of higher-order fluctuations of conserved charges in heavy-ion collisions and the comparison with theory results provided by lattice-QCD simulations and Hadron-Resonance Gas model calculations allow one to get an independent estimate of the chemical freeze-out point [42,43,44], besides the usual one based on the average yields of the various hadronic species. Finally, various combinations of cumulants of conserved charges, allow one to extract information on the nature of charge carriers, and hence of the active degrees of freedom (hadrons or quarks in the case of QCD), in the medium: for a comprehensive review see Ref. [41].
Being the cumulants extensive quantities it is often convenient to take ratios of the latter, so to cancel the dependence on the volume V in Eq. (24), which can be poorly known: in particular, in the case of relativistic heavy-ion collisions, the only indirect information on the size of the medium comes from the multiplicity of produced particles and from HBT correlations. Furthermore one prefers dealing with dimensionless quantities, in which the trivial T 4 dependence of the pressure of a relativistic plasma is factorized. Ratios of cumulants (here, simply of the baryon-number distribution) are then defined in terms of ratios of the following generalized susceptibilities: Here, as independently done in the literature (see e.g. [45]), we focus on the following two interesting ratios of cumulants, directly related to the skewness and the kurtosis of the baryon-number distribution: For the sake of simplicity in the following, with abuse of language, we often refer to the above ratios as normalized skewness and kurtosis. Their interest is first of all related to their behaviour around the chiral transition. The second order cumulant χ B 2 displays a ridge structure in the µ B − T plane along the crossover line [41]. Hence, we expect that χ B 3 changes sign and χ B 4 stays negative around the location of the chiral crossover. This is what one actually observes in Figs Table 4.   Fig. 1. The evolution always starts at time τ0 = 1 fm/c, corresponding to the entropy densities s0 given in Table 4. Fig. 11. The same as in Fig. 10, but for the PQM model.
one varies the center-of-masss energy of the collision [18,19], although no definite conclusions can be drawn so far. As already mentioned, ratios of higher-order cumulants allows one to get information on the nature of the active degrees of freedom in the medium under the different conditions of temperature and density. Consider first the classical limit, in which for the (quasi-)particles of the medium the condition E p − µ T holds. In this case From Eq. (24) it follows that particles obey a Poissonian distributions, with all cumulants equal to the mean, N n c = M . Considering the fluctuations of conserved charges in a relativistic gas one is actually interested in the distribution of the net number of particles (particles minus antiparticles). One gets then the difference of two Poissonian, i.e. a Skellam distribution, for which one has: where the factor (−1) n comes from the multiple derivatives with respect to the chemical potential of the antiparticles, entering with a minus sign. Consider now the following ratio of higher-order cumulants of the net baryon number: where B is the baryon charge carried by the elementary active degrees of freedom. We expect this ratio to be 1 in the confined, chirally-broken phase, where baryon number is carried by particles -the baryons -with B = 1. In Figs. 8-11 we show that in the low-temperature/density phase χ B 3 /χ B 1 and χ B 4 /χ 2 2 actually approach 1 both in the PNJL and PQM models. From Eqs. (17) and (18) one can see that this arises from the role of the Polyakov field, which implements in an effective way colour-confinement suppressing the contribution to the pressure from single-quark (B = 1/3) and two-quark (B = 2/3) clusters, leaving unquenched only the one from states with three quarks of different colours and total baryon number B = 1. We stress that, if the effective valence quarks, with B = 1/3, were active degrees of freedom one would get a much lower value for this ratio, close to 1/9: this is what actually happens in the standard NJL and QM models, which do not implement quark confinement. Notice however that, both in the PNJL and PQM models, in the high-temperature/density limit in which baryon number is carried by quarks, χ B 3 /χ B 1 and χ B 4 /χ 2 2 saturate to a very low value, even smaller than the classical result 1/9 given by Eq. (31), due to the effect of quantum statistics.

Discussion and perspectives
In this paper we explored the thermodynamics of stronglyinteracting matter through two effective Lagrangians developed in the literature to describe the spontaneous breaking of Z 3 (confinement/deconfinement transition) and chiral symmetry: the PNJL and the PQM models. We performed our study in the mean-field approximation, describing the system as a gas of quarks endowed with effective masses and coupled to the Polyakov fields. Both the effective quark masses and the expectation value of the Polyakov fields are obtained requesting the thermodynamic potential to be stationary under variations of the above quantities. Our aim was to obtain very general qualitative information on the phase-diagram of strong interactions and on the behaviour of matter in the different regions of the µ B −T plane, not accessible by lattice-QCD simulations, limited to the case of vanishing or very small baryon density. Our phenomenological motivation was to provide a theoretical guidance for the rich ongoing and future experimental programs at RHIC, SPS, NICA and FAIR, in which heavy-ion collisions at low centerof-mass energy (will) allow the exploration of the chiral/deconfiment transition in the high-density region of the QCD phase-diagram.
The strongest motivation of the above experimental program is the search for the possible critical endpoint in the QCD phase-diagram, were the sharper and sharper crossover would turn into a first-order transition. Hence, we started our analysis identifying the crossover/first-order transition line in the µ B − T plane and establishing the location of the CEP: for (µ CEP B , T CEP ) we obtained the values (875, 121) MeV and (903, 118) MeV in the PNJL and PQM models, respectively. It is of interest to check how close/far these values are from the region currently explored in Beam-Energy Scan ongoing at RHIC. QCD conserves baryon number; furthermore, if dissipative effects are small, entropy production is negligible during the expansion of the fireball arising from the heavy-ion collisions. Hence, the evolution of the medium produced in the collision occurs along trajectories of constant entropy-perbaryon. We estimated the S/BB values of relevance for the BES at RHIC starting from the measured yields of charged particles (assuming S ∼ N ch ) and net protons, obtaining values ranging from 17 to 330 at the lowest ( √ s NN = 7.7 GeV) and highest ( √ s NN = 200 GeV) center-of-mass energy respectively. For comparison, in the PNJL and PQM models, the critical isentropes passing through the CEP correspond to values of (s/n B ) crit of 7.02 and 6.16. If the above estimates are realistic one should conclude that at RHIC we have not yet reached conditions close to the CEP of the QCD phase-diagram. Due to its conceptual and possibly phenomenological importance, in our analysis we explored also the first-order transition region, where two different phases can coexist in equilibrium or -if the evolution of the system is rapid enough compared to the bubble nucleation rate -the system can remain for long in a metastable phase. Hence, in our figures, we also plotted the corresponding spinodal curves predicted by the two models, boundaries of the metastability regions.
If, even at the lowest center-of-mass energy explored in the BES at RHIC, the transition might occur far from the first-order region, are there experimental signatures sensitive to the steeper and steeper crossover as the baryon density increases? We started our investigation considering the speed of sound, whose evolution was studied along the different isentropes. The qualitative behaviour of all the curves referring to the crossover case looks very similar, however -at the quantitative level -at any given time the value of c 2 s in collisions at lower center-of-mass energies is always much lower than the one at higher √ s NN .
Since c 2 s maps energy-density gradients into pressure gradients, responsible for the fluid acceleration, such a sizeable softening of the EoS at lower values of √ s NN should strongly affect the flow of the medium. Experimental signatures of the softening of the EoS can be looked for in the transverse-momentum distributions of the produced hadrons. The average transverse mass m ⊥ of pions, kaons and protons as a function of √ s NN was studied in nucleusnucleus collisions at AGS, SPS, RHIC and LHC. One found m ⊥ − m ≈ 0.2 GeV at the low AGS energies [46], a flat plateau around a higher value at SPS [12] and at the lowest energies of the BES at RHIC [14] and an increasing trend starting from √ s NN ≈ 60 GeV up to LHC energies.
These findings look compatible with the combined effect of the softening of the EoS and of the milder energy-density gradients as √ s NN gets lower. Actually, people suggested that the flat plateau of m ⊥ at SPS energies might be due to the coexistence of the hadronic and deconfined phases typical of a first-order transition. However, more solid conclusions can only be drawn solving the full set of hydrodynamic equations for finite baryon density and with a realistic EoS, with a steeper and steeper crossover eventually turning into a first-order transition.
If the speed of sound mainly affects the momentum distribution of the produced hadrons, the fluctuations of their yields -more precisely the ones associated to conserved quantities, like net protons, kaons and charged particlesare related to the higher-order susceptibilities of baryonnumber, strangeness and electric charge. Accordingly, we addressed their evaluation starting from the mean-field thermodynamic potential of the PNJL and PQM models, focusing in this paper only on baryon-number fluctuations. At the CEP the order parameter (the σ field or chiral condensate) is characterized by an infinite correlation length ξ leading to a divergence of its cumulants and of all the cumulants of quantities coupled to the latter, like the baryon number. Notice that in heavy-ion collisions, due to the finite size and lifetime of the medium, one would not observe any actual divergence of the correlation length and it is thus necessary to focus on higherorder cumulants, which display a stronger sensitivity on ξ. Hence, we studied the ratios of the generalized baryonnumber susceptibilities χ B 3 /χ B 1 and χ B 4 /χ B 2 along lines of fixed s/n B of interest for the BES at RHIC, finding that for the trajectories passing closer to the CEP -for which the crossover is steeper -the above quantities display sizeable oscillations and rapid changes of sign in the transition region. Clearly, one aims at finding signatures of the above non-trivial behaviour in the proximity of the CEP in the experimental fluctuations of conserved charges (baryon number, strangeness and electric charge), looking for deviations from a trivial Skellam statistics, although how far the chemical freeze-out point lies from the chiral transition plays a crucial role. For the moment the above higher-order cumulants were used for a less ambitious purpose, namely to get an independent estimate of the chemical freeze-out point in the µ B −T plane, besides the usual one based on the average yields of identified hadrons.
Our work has to be considered just as a first step towards the aim of providing solid theoretical qualitative guidance within chiral effective models for the experimental exploration of the phase-diagram of stronglyinteracting matter. We plan to improve our models in various aspects. First of all we plan to go beyond the present mean-filed approximation, largely employed in the literature, but unable to provide a realistic picture of the hadronic phase, where pions, kaons, protons... are important degrees of freedom contributing to the thermodynamics. Notice that mesons (and also baryons) can be obtained in the PNJL model as solutions of a Bethe-Salpeter equation [47,48,49], but one should go one step further and add self-consistently their contribution to the thermodynamic potential [50]. Secondly, we plan to add to the models a vector interaction. We should get then a more realistic description of the phase diagram, where presently the value of the nuclear-matter density n B ≈ 0.16 fm −3 lies in the phase-coexistence region, representing an unphysical feature of the two models. Finally, we plan to carry out a deeper study of the first-order transition, investigating aspects like the interface energy between the two phases, the rate of bubble/droplet nucleation and the possible occurrence of spinodal instabilities under the experimental conditions. We leave the above items for forthcoming publications.