The presence of a critical end point in the QCD phase diagram from the net-baryon number fluctuations

The net-baryon number fluctuations for three-flavor quark matter are computed within the Polyakov extended Nambu$-$Jona-Lasinio model. Two models with vanishing and nonvanishing vector interactions are considered. While the former predicts a critical end point (CEP) in the phase diagram, the latter predicts no CEP. We show that the nonmonotonic behavior of the susceptibilities in the phase diagram is still present even in the absence of a CEP. Therefore, from the nonmonotonic behavior of the susceptibilities solely, one cannot assume the existence of a CEP. We analyze other possible properties that may distinguish the two scenarios, and determine the behavior of the net-baryon number fluctuations and the velocity of sound along several isentropes, with moderate and small values. It is shown that the value of the susceptibilities ratios and the velocity of sound at two or three isentropic lines could possibly allow to distinguish both scenarios, a phase diagram with or without CEP. Smoother behaviors of these quantities may indicate the nonexistence of a CEP. We also discuss the critical behavior of the strange sector.


I. INTRODUCTION
The quest for the QCD chiral critical end point (CEP) in the phase diagram, together with the nature of the phase transition between hadron matter and quark gluon plasma (QGP), are open questions that have attracted the attention of the physical community for some years [1]. Remarkable theoretical and experimental efforts [2] are being made to unveil the rich details of the QCD phase structures [3]. Experimentally, one of the main goals of the heavy ion collision (HIC) program is the possible existence and location of the CEP on the QCD phase diagram, with great developments over the last years [4][5][6][7].
In relativistic HIC, the measurement of fluctuations of conserved quantities, such as baryon, electric charge, and strangeness number, play a major role in the experimental search for the CEP. Indeed, experimental measurements of cumulants of net-proton (proxy for net-baryon), net-charge, and net-kaon (proxy for net-strangeness) are expected to carry significant amounts of information on the medium created by the collision (for a review, see [8][9][10][11]). Fluctuations are studied by measuring event-byevent fluctuations: a given observable is measured on an event-by-event basis and its fluctuations are studied for the ensemble of events [10].
Particularly relevant are the cumulants of the netbaryon number because a second-order phase transition occurs at the CEP, resulting in divergences of correlation lengths for a static system of infinite size. The cumulants of the baryon number thus diverge at the CEP [12,13]. The study of the kurtosis [14] and the skewness [15] for the net-baryon number fluctuation distributions is essential as they are related to higher-order cumulants that can be extracted from event-by-event fluctuations in HIC experiments. Once they are constituted by ratios of cumulants they are independent of the volume of the system.
Other models have been employed to study higherorder baryon number susceptibilities at finite temperature and density like the Polyakov-loop extended quarkmeson model [25], where the influence of repulsive vectorinteractions on this fluctuations was also analyzed, the hybrid quark-meson-nucleon model [26], or the SU(3) flavor parity-doublet quark-hadron model [27] where the higher-order baryon number susceptibilities near the chiral and the nuclear liquid-gas transitions were investigated.
The eventual location of the CEP can be affected by several conditions such as the presence of external magnetic fields or the strangeness and isospin content of the medium [28][29][30]. The study of the CEP location has been undertaken using different versions of the NJL and PNJL models. In particular, it was shown that the presence of repulsive vector interactions affects strongly the position of the CEP. The role played by them were analyzed in detail in [29,31]. The calibration of these models at high densities requires the existence of experimental data or neutron star observables. Particularly relevant is the introduction of repulsive interactions, namely the vectorisoscalar terms, that seems to be necessary to describe 2M hybrid stars [32].
The chiral restoration of strange quarks may play an important role inside neutron stars. In particular, if this transition occurs at densities that can be found inside compact stars, pure quark matter [32], or, exotic quark phases such as the color-flavor-locked (CFL) phase could be realized in their interior [33]. Besides, a phase transition could also have an important effect on the mean-free path of neutrinos in a protoneutron star as discussed in [34]. The cooling of protoneutron stars during the first seconds is essentially driven by the neutrinos that diffuse out of the star. A phase transition would give rise to a opalescence like phenomena reducing a lot the neutrino mean-free path, and, therefore, allowing for a much larger interaction of neutrinos with matter.
In this work, we study the phase diagram using the (2+1)-flavor Nambu-Jona-Lasinio model coupled to the Polyakov loop, designated as PNJL model, from the point of view of the kurtosis and skewness of net-baryon number fluctuations. It is expected that in HIC the fireball evolves along isentropes, lines with constant entropy per baryon, and, therefore, we analyze how these quantities, as well as the velocity of sound, behave along isentropes. Our main objective is to identify the similarities and differences of a QCD phase diagram which has a CEP or not, namely when a sufficiently strong repulsive vector interaction is taken into account.
The model is succinctly reviewed in Sec. II, while the results are discussed in Sec. III. Finally we draw our conclusions in Sec. IV.

II. MODEL AND FORMALISM
The Lagrangian density for the Polyakov extended Nambu-Jona-Lasinio (PNJL) model reads where the quark field is represented by q = (u, d, s) T in flavor space, andm c = diag f (m u , m d , m s ) is the corresponding (current) mass matrix. The L sym and L det denote, respectively, the scalar-pseudoscalar and the 't Hooft six-fermion interactions [35,36], The vector interaction is given by [37] The effective gluon field is given by A µ = g strong A µ a λa 2 , where A µ a represents the SU c (3) gauge field. The spatial components are neglected in Polyakov gauge at finite temperature, i.e., A µ = δ µ The Polyakov loop value is defined as the trace of the Polyakov line, which is the order parameter of the Z 3 symmetric/broken phase transition in pure gauge. For the pure gauge sector we use the following effective potential [38], Its parametrization values are a 0 = 3.51, a 1 = −2.47, a 2 = 15.2, and b 3 = −1.75 [38], while the critical temperature is set to T 0 = 210 MeV. The divergent ultraviolet sea quark integrals are regularized by a sharp cutoff Λ in three-momentum space. For the NJL model parametrization, we consider: Λ = 602.3 MeV, m u = m d = 5.5 MeV, m s = 140.7 MeV, G s Λ 2 = 1.835, and KΛ 5 = 12.36 [39].
Fluctuations of conserved charges, such as the baryon number, provide vital information on the effective degrees of freedom and on critical phenomena. They behave characteristically in a thermal equilibrium medium. If there is a CEP in the phase diagram of strongly interacting matter, these fluctuations are then expected to provide characteristic signatures that, hopefully, can be experimentally observed. For a static system of infinite size, the fluctuations of baryon number diverge at the CEP (second-order phase transition point). However, the created medium in HIC experiments has both finite size and lifetime that restricts its correlation length and, instead of divergent fluctuations, only moderate enhancements are expected. Fluctuations of conserved charge are characterized by their cumulants or susceptibilities. The present work focuses on the baryon number charge susceptibilities. The nth-order net-baryon susceptibility is given by Different susceptibility ratios χ n B (T, µ B )/χ m B (T, µ B ) are calculated in order to eliminate the volume dependence, allowing for a comparison with experimental observables.
In this work, we analyze the following ratios where M = V T 3 χ 1 B is the mean, σ 2 = V T 3 χ 2 B the variance, S B the skewness, and κ is the kurtosis of the netbaryon number distribution.

III. RESULTS
We analyze, herein, the net-baryon susceptibilities on the (T, µ B ) plane. Two PNJL models are analyzed: (i) a model with no vector interactions G V = 0, which predicts a CEP; and (ii) a model with vector interactions G V = 0.72G s , which predicts no CEP. We want to discuss what distinguishes these two scenarios. In the following, symmetric quark matter is considered: where µ i are the chemical potential of each quark flavor and µ B is the baryonic chemical potential.
The hydrodynamical expansion of a HIC fireball is expected to follow trajectories of constant entropy per baryon, s/ρ B , known as isentropes. These trajectories contain important information on the adiabatic evolution of the system. It is thus interesting to analyze the susceptibility ratios [Eqs. (7)] along different isentropes [40]. It is important to note that while the net charge and the net strangeness are not constrained in the present work; in a HIC, however, the ratio of electric charge over baryon number is Q/ρ B 0.4 and no net strangeness is produced, n s = 0.
The phase diagrams for the chiral and deconfinement transitions are presented in Fig. 1 while the G V = 0.72G s model has no CEP, and the (approximate) chiral restoration is thus attained via an analytic transition (crossover) over the whole phase diagram. The chiral (dashed line) and deconfinement (dotted line) crossover boundaries are determined by the location (T, µ B ) of the maximum of the order parameter susceptibilities (the point where fluctuations are largest). It is interesting that the crossover boundaries show similar behavior for both models: the gap between the deconfinement and chiral crossovers reduces with increasing µ B and becomes zero for some µ B values, which turns out to be near the CEP for G V = 0, above which they separate and follow distinct paths.
Both boundaries, the chiral phase transition and the deconfinement phase transition boundaries, are determined from the peaks of the susceptibility. The crossing of the deconfinement and chiral phase transitions has already been observed before [41] and it is possible to identify the crossing from the calculation of the susceptibility peaks at fixed temperatures: before and after the crossing they are two distinguishable peaks. At the crossing, that stretches along a finite range of temperatures, the two peaks overlap. The crossing region includes part of the chiral crossover for both models, and, in the case of the model with a CEP, also the CEP, and part of the first-order phase transition.
We thus conclude that, due to the mixing between the gluonic and quarkionic degrees of freedom, the chiral phase transition has a strong influence on deconfinement transition. This is reflected on the behavior of the deconfinement transition at the light quark and the strange quark chiral transition. At the light quark transition, the crossing temperature is not much affected, but the crossing chemical potential is tightly connected with the position of the chiral transition and the crossing follows, as referred above, the chiral crossover or both the chiral crossover and first-order transition. As a consequence, the crossing occurs at a much larger chemical potential for the G V = 0.72G s model. A similar interconnection is observed at the strange chiral crossover in Fig. 2   two last trajectories cross the crossover line above the CEP of the G V = 0 model. The first three values allows us to discuss the phase diagram at low T and high µ B , where the (approximate) chiral restoration of the strange quark occurs. For the G V = 0 model, the susceptibilities exhibit a nonmonotonic dependence near the CEP, whose behavior strongly depends on the direction on which the CEP is approached. The susceptibilities diverge at the CEP, with the divergence being stronger as higher susceptibilities orders are considered. An interesting result is present at low T and high µ B . Despite the transition for the strange quark being just a crossover, and, therefore, without any nonanalytic behavior, a similar CEP structure is seen at µ B ≈ 1500 MeV for the susceptibilities. This indicates that a slight change on the model parametrization might induce a first-order phase transition for the strange quark, and a corresponding CEP. The χ 3 B and χ 4 B values for the G V = 0.72G s model show precisely this behavior for the light quark sector: even though there is no CEP, and the chiral transition occurs via a crossover over the whole phase diagram, the nonmonotonic behavior of the susceptibilities is still present, as discussed within the NJL model [18]. The study of a scenario with a hypothetical negative T CEP for the light CEP, obtained by varying the value of the anomaly-induced six-fermion term, K, was done in [42]; it was shown that the magnitude of the susceptibilities also changes significantly if a hypothetical negative temperature CEP is taken into account.
The ratios χ 4 B /χ 2 B and χ 3 B /χ 1 B are shown in Fig. 3. The sudden decrease near the deconfinement pseudocritical temperature (dotted black line) indicates that both quantities are valuable signatures of deconfinement transition. As noted in [19], the statistical confinement, provided by the Polyakov loop (at low temperatures, when Φ,Φ → 0, contributions coming from one-and two-quark states are suppressed, while three-quark states are not [43]), is essential to obtain a low-temperature limit for the susceptibility ratios that is consistent with the hadron resonance gas model. The results for the G V = 0.72G s model clearly show that the nonmonotonic behavior of χ 3 B and χ 4 B , which signals the presence of a critical behavior, is still present even in the absence of a CEP. The nonmonotonic behavior persists, with a smaller intensity, up to almost the same temperature as for the G V = 0 model. To make this feature clear, we show the negative region of χ 4 B /χ 2 B in Fig. 4. Despite the strong vector interaction used, it is remarkable that, for the G V = 0.72G s model (red), the χ 4 B /χ 2 B < 0 region extends from zero up to temperatures similar with the ones obtained for G V = 0 (blue). This indicates that, at higher temperatures, both models are not discernible exclusively from the sign change of χ 4 B /χ 2 B . Actually, a region with χ 4 B /χ 2 B < 0 is still present (at lower temperatures though) even when the vector interaction strength is increased up to G V ≈ 1.4G s . If instead of looking at the whole negative region of χ 4 B /χ 2 B < 0, one considers the stronger fluctuation region χ 4 B /χ 2 B < −200 the following pattern is seen: while this region extends to a range of ∆µ B ≈ 100 MeV and ∆T ≈ 20 MeV for G V = 0, we get ∆µ B ≈ 20 MeV and ∆T ≈ 60 MeV for G V = 0.72G s . These different ranges on T and µ B for the two different scenarios could help distinguish them, taking only into account the behavior of the fluctuations. It should, however, be recalled that if the CEP exists, for moderate temperatures and high enough baryonic density the line of first-order transition could be crossed during the evolution of the fireball, giving rise to effects like multifragmentation [10,44]. This would be a region where our no CEP model would present fluctuations similar to the ones existing in a model with CEP, above the CEP. To complete the discussion, in the following we analyze the isentropic trajectories with a small s/ρ B within the two scenarios.
The comparison of the isentropic trajectories between G V = 0 (solid lines) and G V = 0.72G s (dashed lines) models is in Fig. 5. The trajectories differ for high values of T and µ B , i.e., as soon as the system becomes denser enough for the vector interactions to set in. Two features that distinguish the G V = 0 model from the G V = 0.72G s is the behavior of the trajectories near the CEP and the existence of a unstable spinodal region. The trajectories with low s/ρ B values get enclosed into the unstable spinodal region when crossing the first-order phase transition to the chiral broken phase. As the system enters into the unstable spinodal region, the rapid formation of fragments of high density matter that occur should enhance the baryon number fluctuations [10]. Due to the absence of spinodal region for the G V = 0.72G s model, such effect does not occur and the susceptibilities have an analytic behavior.
In Fig. 6, we show the χ 4 B /χ 2 B (top) and χ 3 B /χ 1 B (bottom) values along the s/ρ B = 14 (left) and the s/ρ B = 10 (right) isentropes (these isentropic trajectories are shown in Fig. 2 and 3). As the value s/ρ B of the isentropic trajectory decreases, we are covering a higher µ B region on the phase diagram. As we move from s/ρ B = 14 to s/ρ B = 10, we are then approaching a region of higher baryon fluctuations that reflects the vicinity of a CEP. While the fluctuations of χ 4 B /χ 2 B and χ 3 B /χ 1 B grow with decreasing s/ρ B , they also become constrained to a smaller temperature region (this is clear through the shape of the blue region in Fig. 4). The decreasing gap between the chiral and deconfinement transitions with increasing µ B , which vanishes at the CEP (see Figs. 2 and 3), is also reflected in the fluctuations: for s/ρ B ≥ 14 a two peak structure is present on the left side of the fluctuation (for s/ρ B = 14, a small bump at T ≈ 150 MeV is barely seen).
This two peak structure, which reflects the deconfinement/chiral restoration gap, is clearer when a vector interaction is included. The fluctuations for the G V = 0.72G s model, over the same isentropic trajectories, are shown in Fig. 7. The fluctuations along the s/ρ B = 14 trajectory show a two peak structure again on the left side of the fluctuation. Despite the existence of a sign change of χ 4 B /χ 2 B (top) for both models (also seen in Fig. 4), their intensity is weaker for the G V = 0.72G s model, allowing one to notice the effect of the deconfinement transition on the fluctuations ratios.
Let us now focus on the crossover region at low temperatures, i.e., low s/ρ B values, for the G V = 0.72G s model. In Fig. 8  particles originating from the fluid elements at the freezeout stage [45]. The values of square sound velocity are extracted from the widths of rapidity distributions [45][46][47]. For example, from the measured data on the widths of the pion rapidity spectra, v 2 s in the dense stage of the reactions has been extracted [48].
In Fig. 9, we show v 2 s along s/ρ B = 10 (top) and s/ρ B = 1 (bottom) for G V = 0 (blue) and G V = 0.72G s (red) models. As the isentropic trajectories follow specific paths, (T, µ B ), on the phase diagram (see Figs. 2 and 3), we show the v 2 s dependence on temperature (left) and baryon chemical potential (right), for both isentropic trajectories. For each isentropic, we give the values of v 2 s (T, µ B ) at its minimum and at the chiral and deconfinement boundaries in Table I.
Considering in first place the s/ρ B = 10 trajectory as a function of T (Fig. 9, left upper panel), we see that the minimum of v 2 s (T, µ B ) occurs closer to the deconfinement pseudocritical boundary (circles) than to the chiral pseudocritical boundary (squares). While the temperature dependence of the s/ρ B = 10 isentrope is a single-valued function, the same does not hold for its µ B dependence ( Fig. 9, right upper panel). The loop behavior for the G V = 0 model (blue curve in upper right panel of Fig. 9) rises from the bending effect towards the CEP that the s/ρ B = 10 isentrope undergoes when crossing into the chiral broken region (solid red line in Fig. 5). This effect occurring in v 2 s can then be seen as a signal for the vicinity of a CEP (if some kind of bending effect into the CEP exists), once this effect is not seen for the G V = 0.72G s model. For s/ρ B = 1 (lower panels of Fig. 9), the v 2 s shows negative values for G V = 0 model (blue curve), reflecting the first-order phase transition that occurs at lower T . It is interesting to note that, for small values of s/ρ B , the local minimum of v 2 s at µ B ≈ 1480 MeV is associated with the crossover of the strange quark.

IV. CONCLUSIONS
We have analyzed the net-baryon number fluctuations for three-flavor quark matter within the Polyakov extended Nambu-Jona-Lasinio model. For a strong enough vector interaction intensity, the model predicts no CEP in the phase diagram. From the net-baryon number fluctuations one concludes that, even in the absence of a CEP, the nonmonotonic behavior persists. Therefore, the existence of a CEP cannot be taken solely from the existence of nonmonotonic behavior on the net-baryon number susceptibilities.
We have analyzed further other possible properties that may distinguish the two scenarios: for the no CEP model (G V = 0.72G s ), large fluctuations in the susceptibility ratios occur only at small T , and the values of v 2 s are almost unchanged at moderates s/ρ B values. The values of the susceptibility ratios along two or three isentropic lines would possibly allow us to distinguish both cases. Also, the value of the sound velocity at the chiral transition for two or three isentropes would give some useful information. For the G V = 0 model, by going from s/ρ B = 9 to 14, the value of v 2 s increases at least 50% for each step ∆(s/ρ B ) = 2. Instead, for the no CEP model, the change is of the order of 10%. It should be noticed, however, that in the present work we have discussed infinite size matter. For a finite system it is expected that the signals we have discussed are less intense but still might allow to distinguish both scenarios.
We have shown that, for high chemical potentials and low temperatures, a signature of the strange quark chiral symmetry restoration is observed in a decrease of the sound velocity and in a region with negative χ 4 B /χ 2 B . Acknowledgments the CENTRO2020 program. Partial support comes from "THOR" (COST Action CA15213) and "PHAROS" (COST Action CA16214).
which satisfy the following gap equations: The values of Φ andΦ are the solutions of The thermodynamical quantities are determined via the thermodynamical potential (see [40]). The pressure is given P (T, µ i ) = −Ω(T, µ i ), the density of the i quark, ρ i , is given by while the the entropy density, s, is given by The energy density, E, comes from the following fundamental relation of thermodynamics