The kurtosis of net baryon number fluctuations from a realistic Polyakov–Nambu–Jona-Lasinio model along the experimental freeze-out line

Firstly we qualitatively analyze the formation of the dip and peak structures of the kurtosis κσ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa \sigma ^2$$\end{document} of net baryon number fluctuation along imagined freeze-out lines and discuss the signature of the existence of the QCD critical end point (CEP) in the Nambu–Jona-Lasinio (NJL) model, Polyakov-NJL (PNJL) model as well as μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document}-dependent PNJL(μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} PNJL) model with different parameter sets, and then we apply a realistic PNJL model with parameters fixed by lattice data at zero chemical potential, and quantitatively investigate its κσ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa \sigma ^2$$\end{document} along the real freeze-out line extracted from experiments. The important contribution from gluodynamics to the baryon number fluctuations is discussed. The peak structure of κσ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa \sigma ^2$$\end{document} along the freeze-out line is solely determined by the existence of the CEP mountain and can be used as a clean signature for the existence of CEP. The formation of the dip structure is sensitive to the relation between the freeze-out line and the phase boundary, and the freeze-out line starts from the back-ridge of the phase boundary is required. To our surprise, the kurtosis κσ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa \sigma ^2$$\end{document} produced from the realistic PNJL model along the experimental freeze-out line agrees with BES-I data well, which indicates that equilibrium result can explain the experimental data. It is worth to point out that the extracted freeze-out temperatures from beam energy scan measurement are indeed higher than the critical temperatures at small chemical potentials, which supports our qualitative analysis.


Introduction
The phase transition and phase structure of Quantum Chromodynamics (QCD) under extreme conditions is the main topic of relativistic heavy ion collisions, and it is also highly related to the evolution of the early universe and the equation of state inside the compact stars. Lattice QCD calculation shows that at small baryon density, the QCD phase transitions including the chiral phase transition as well as deconfinement phase transition are of crossover at finite temperature [1][2][3]. From symmetry analysis and effective chiral models, it is generally believed that at high baryon density the chiral phase transition is of first order and there exists a QCD critical end point (CEP) for chiral phase transition in the temperature and baryon chemical potential plane [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. The QCD CEP have been widely analyzed in different models, e.g., Nambu-Jona-Lasinio (NJL) model, the Polyakov-loop improved NJL (PNJL) model, linear sigma model, quarkmeson (QM) model, the Polyakov-loop improved QM model, the Dyson-Schwinger equations (DSE), and the holographic QCD model [6][7][8][9][10][11][12][13][14][15][16][17][18][19]24]. However, different models even the same model with different parameter sets give various location of CEP [25]. Therefore, to search for the existence of the CEP and further to locate the CEP is one of the most central goals at Relativistic Heavy Ion Collisions (RHIC) as well as for the future accelerator facilities at Facility for Antiproton and Ion Research (FAIR) in Darmstadt and Nuclotron-based Ion Collider Facility (NICA) in Dubna.
The cumulants of conserved quantities up to fourth order of net-proton, net-charge and net-kaon multiplicity distributions have been measured in the first phase of beam energy scan program (BES-I) at RHIC for Au+Au collisions at √ s N N = 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and 200 GeV, and the results are summarized in [26][27][28]. A non-monotonic energy dependent behavior for the kurtosis of the net pro-ton number distributions κσ 2 has been observed in the most central Au+Au collisions: κσ 2 firstly decreases from around 1 at the colliding energy √ s N N = 200 GeV to 0.1 at √ s N N = 20 GeV and then rises quickly up to around 3.5 at √ s N N = 7 GeV. From the BES-I observed oscillation behavior of kurtosis, we need to know what information we can extract about the CEP and where is the CEP located. Before the second phase of beam energy scan (BES-II) at RHIC performed in 2019-2020, it is very urgent for both experimentalists and theorists to extract a clean signature to identify the existence of the QCD CEP and further to locate the CEP. On the one hand, we should try to extract useful information about QCD phase transitions from the measurement along the freeze out line, and to find the evidence of the QCD CEP. On the other hand, we should explore carefully how QCD phase transitions will shed light on properties of the cumulants of conserved quantities at freeze-out.
By using the Ising model, it was pointed out in [21] that the quartic cumulant (or kurtosis) of the order parameter is universally negative when the critical point is approaching to the crossover side of the phase separation line. From the results in [21], the oscillation behavior of the kurtosis along the freeze out line has been regarded as a typical signature of the existence of the CEP. Many interests have focused on the sign changing of various cumulants around the CEP, and the sign changing for higher order susceptibilities has been recently discussed in [11]. It is noticed that the sign changing for the 6th and 8th order susceptibilities starts at the baryon chemical potential quite far away from the CEP, it may indicate that the sign changing of cumulans of conserved quantities is not directly related to the CEP. In the holographic QCD model [29], we explained that the sign changing of the baryon number susceptibilities along the freeze-out line is not necessarily related to the CEP, but the peaked baryon number susceptibilities along the freeze-out line is solely determined by the CEP thus can be used as an evident signature for the existence of the CEP, and the peak position is close to the location of the CEP in the QCD phase diagram. Furthermore, it is found that at zero chemical potential, the magnitude of κσ 2 around the phase transition line in the gluodynamics dominant holographic QCD model is around 1, which is in agreement with lattice result, and is much larger than that in the Ising model (close to zero), and also much larger than that in the NJL model (around 0.1). Therefore, it is natural to speculate that the gluodynamics contribution is dominant to the baryon number susceptibilities, which is quite surprising! In order to check the gluodynamics contribution to the baryon number susceptibilities, in this work, we will try to explore the structure of the kurtosis of the net proton number distribution κσ 2 along the freeze out line in the Nambu-Jona-Lasinio (NJL) model, the Polyakov-loop improved NJL (PNJL) model as well as μ-dependent Polyakov-loop potential improved NJL (μPNJL) model. Also in order to investi-gate the relation between the CEP and the structure of the κσ 2 along the freeze out line, we will need to shift the location of the CEP by introducing the interaction in the vector channel. This paper is organized as following: After Introduction, in Sect. 2, we give a brief introduction to the two-flavor NJL model, the Polyakov-loop improved NJL (PNJL) model as well as μ Polyakov-loop improved NJL (μPNJL) model, and qualitatively analyze the formation of the dip and peak structures of the kurtosis κσ 2 of net baryon number fluctuation along imagined freeze-out lines and discuss the signature of the existence of the QCD critical end point (CEP). Then in Sect. 3 we apply a realistic PNJL model with parameters fixed by lattice data at zero chemical potential, and quantitatively investigate its κσ 2 along the real freeze-out line extracted from experiments. Finally, the discussion and conclusion part is given in Sect. 4.

Qualitative analysis of baryon number fluctuations in the NJL, PNJL and μPNJL models
In order to qualitatively analyze the formation of the dip and peak structure of the kurtosis of the net baryon number distribution κσ 2 along the freeze-out line, as well as to investigate the gluodynamics contribution to κσ 2 , in this section we will compare κσ 2 in the framework of NJL model, the Polyakov-loop improved NJL (PNJL) model as well as NJL with μ-dependent Polyakov-loop potential (μPNJL) model. For each model, except the coupling constant in the vector channel, the parameters from quark part are fitted by vacuum properties as well as pion mass and decay constant, and the parameters from the Polykov-loop potential part are fitted from the Lattice results of the equation of state at μ = 0. We will intendedly shift the location of the CEP in the models by changing the coupling constant in the vector channel, to check the peak structure of κσ 2 along the imagined freezeout lines and its relation with the location of CEP.

The NJL, PNJL and μPNJL models with vector interaction
In order to shift the location of the CEP, we introduce the two-flavor NJL model with the vector interaction, and the Lagrangian is given by [30][31][32] where ψ = (u, d) T is the doublet of the two light quark flavors u and d with the current mass m = m u = m d , and τ = (τ 1 , τ 2 , τ 3 ) the isospin Pauli matrix. G S and G V are the coupling constants in the (pseudo)scalar channel and the vector channel, respectively. By introducing the auxiliary fields for scalars and vectors, and in the vacuum we take the meanfield approximation with the quark condensate and the net quark number density of flavor i respectively. Then the thermodynamical potential of this NJL model takes the following form: with N c = 3 the number of colors. The quark quasiparticle energies E i and constituent quark masses M i for flavors i = (u, d) are given by and the effective chemical potentials are shifted bỹ In order to solve the minimum of the thermal potential NJL , we have the following gap equations and For numerical calculations, we fix m u = m d = 5.5 MeV, N c = 3, N f = 2 in all the models. In the regular NJL model, we choose the parameters as in Ref. [40] by fitting the pion mass and decay constant. To shift the location of the CEP in the NJL model we choose different coupling constant in the vector channel: G V = −0.5G S , G V = 0 and G V = 0.67G S [11], which is indicated as NJL-1, NJL-2, NJL-3 respectively, and three sets of parameters are shown in Table  1 In the 2-flavor PNJL and μPNJL models, we need to add the Polyakov-loop effective potential U( ,¯ , T ) with the following ansatz [33][34][35][36] U( ,¯ , T ) where 3 . The critical temperature of the confinement-deconfinement phase transition T 0 in pure gluon system is 270 MeV, and will be rescaled to about 220 MeV because of the presence of fermion fields. The thermodynamical potential of this PNJL model is given by with the same definitions of E i , M i andμ i as in the NJL model. The gap equations are determined by and In the PNJL model, we fix T 0 = 270 MeV, and choose the parameter sets by fitting the experimental values of pion decay constant f π = 92.3 MeV and the pion mass m π = 139.3 MeV when G V = 0. And we choose different G V for PNJL-1, PNJL-2 and PNJL-3 as shown in Table 2. For the PNJL-1, at zero chemical potential μ = 0, the critical temperature for the chiral phase transition is T      Table 3 Three parameters sets used in the μPNJL model, and the corresponding critical temperatures and chemical potentials at the critical end point The thermodynamical potential of μPNJL model μPNJL has the same form of PNJL model, and the gap equations take the following form: and In the μPNJL model, we also choose three different parameter sets as shown in Table 3 [39]. For the μPNJL-1, at μ = 0, the critical temperature for the chiral phase transition is T with the pressure P = − which is just the minus thermodynamical potential. The cumulants of baryon number distributions are given by By introducing the variance σ 2 = C B 2 , kurtosis κ = C B 4 (σ 2 ) 2 , one can have the following relation between observable quantities and theoretical calculations: which relates κσ 2 with the ratio of fourth and second order cumulants of net-baryon number fluctuations. Because the measurement of κσ 2 from BES-I is along the freeze-out line, here we define three different imagined freeze-out lines for the NJL, PNJL and μPNJL models as: Here I = 1, 2, 3 is the label of different model e.g. NJL (PNJL or μPNJL)-I and p, μp is for the PNJL and μPNJL Fig. 1 The kurtosis of baryon number fluctuation κσ 2 as a function of the temperature at zero baryon chemical potential in the NJL model, PNJL model and μPNJL model, and compare with the lattice result in [41]. κσ 2 = 1 in the hadron resonance gas (HRG) limit and κσ 2 0.068 in the ideal free quark gas (FQG) limit model, respectively. And T I p (μ B ) is the function fitting of chiral phase transition line and both T and μ B are in MeVs. These three imagined freeze-out lines just qualitatively indicate how far the freeze-out line is away from the phase boundary, which is an important information to analyze the structure of κσ 2 along the freeze-out line.
In Fig. 1, we show κσ 2 as a function of T /T χ 0 in three parameter sets of NJL model, and three parameter sets of PNJL model and three parameter sets of μPNJL model at zero chemical potential, and compare with lattice results in [41], with T χ 0 the critical temperature for chiral phase transition at zero chemical potential μ B = 0. The critical temperature of lattice data at μ B = 0 is 154 ± 9 MeV and we choose 154 MeV. The value of κσ 2 of net baryon number fluctuations is unity in the limit of hadron resonance gas (HRG) and, in the ideal free quark gas (FQG) limit at infinite temperature it takes the value of κσ 2 0.068, which are also shown in Fig. 1. It is worthy of mentioning that from the lattice result shown in Ref. [41], that at the critical temperature T χ 0 , the magnitude of κσ 2 is around 0.8, which is smaller than its value of 1 at the temperature of T dec 0 140 MeV, at which quark matter transfers to hadron gas. Here T dec 0 describes a physical "confinementdeconfinement" transition, which is different from the T D 0 describing the confinement-deconfinement phase transition by using the order parameter of Polyakov loop.
In the NJL model, the magnitude of κσ 2 at T χ 0 is around 0.4, 0.2, 0.1 in the NJL-1, NJL-2 and NJL-3, which is much smaller than the lattice result. In the whole temperature region T > T 0 , the magnitude of κσ 2 is quite small comparing with the lattice results.
When the gluodynamics contribution is taken into account in the PNJL model and μPNJL model, the critical temperatures for the chiral phase transition T χ 0 and deconfinement phase transition T D 0 are separated, and the critical temperature for the deconfinement phase transition T D 0 is lower than that of the chiral phase transition T χ 0 at zero chemical potential. It is observed that the magnitude of κσ 2 shows a peak at the critical temperature for deconfinement phase transition T D 0 , and in the hadron gas phase with T < T D 0 , the magnitude of κσ 2 decreases with the decreasing temperature, and matches the HRG limit at low temperature.
From the magnitude of κσ 2 in the NJL model, PNJL model and μPNJL model, we can see that the dominant contribution to κσ 2 comes from gluodynamics at zero chemical potential. In the future for model construction, κσ 2 as a function of the temperature can be used to constrain models.

The kurtosis of baryon number fluctuation κσ 2 in the NJL model
In Fig. 2, we show the 3D plot for the kurtosis of baryon number fluctuation κσ 2 as a function of the temperature and baryon chemical potential in the NJL model. In Fig. 3 we show the chiral phase transition line and 2D plot for κσ 2 as a function of the baryon chemical potential along different freeze-out lines for NJL-1,NJL-2 and NJL-3, respectively. By changing the coupling constant in the vector channel, the location of the CEP will shift. The CEP in the NJL-1 and NJL-2 is located at (μ E B = 796.7 MeV, T E = 76.8 MeV), (μ E B = 1005.2 MeV, T E = 34.6 MeV), respectively, and it is observed that κσ 2 develops a high CEP mountain in the NJL-1 and NJL-2. There is no CEP in the NJL-3, therefore, no CEP mountain of κσ 2 develops in the NJL-3. From the 3D plot in Fig. 2, one can observe an obvious chiral phase boundary in the crossover side, and the magnitude of its ridge decreases with the increase of the baryon chemical potential.
Here we define the ridge as along the phase boundary, the back/front ridge as the higher/lower temperature side comparing with the phase boundary, respectively. If there is a CEP located on the phase diagram, a CEP mountain rises up around the CEP, and one can observe a negative region of κσ 2 around the CEP above the chiral phase boundary extended from the crossover side.
Because the measurement of heavy-ion collision is along the freeze-out line, we choose three different freeze-out lines It is observed that when the freeze-out line crosses the CEP mountain or crosses the foot of the CEP mountain, there will be a peak showing up for κσ 2 along the freeze-out line, and the location of the peak is close to the location of the CEP mountain. In the case of no CEP, it is found that κσ 2 keeps flat in almost the whole region and does not show any structure. In the NJL-1 and NJL-2 with the existence of CEP in the phase diagram, it is found that for the first case freeze-out line, because the freeze-out line starts from the back ridge of the chiral phase boundary and goes through the negative region of κσ 2 , then crosses the CEP mountain, therefore, κσ 2 decreases from around 0.2 and then down to negative value then rises up quickly and shows up a high peak around the critical chemical potential, which shows a dip and then a peak structure of κσ 2 along the freeze-out line. For the second case of freeze-out line, because the freeze-out line starts from the chiral phase boundary, and then crosses the foot of the CEP mountain, we can only see a peak structure of κσ 2 along the freeze-out line at high chemical potential. For the third case of freeze-out, if the freeze-out line is far away from the phase boundary, one can only observe a weak peak of κσ 2 along the freeze-out line.
It is found that the magnitude of κσ 2 in the NJL model at small baryon chemical potential region is small comparing with experiment measurement which is around 1.

The kurtosis of the baryon number fluctuation κσ 2 in the PNJL model
In Fig. 4, we show the 3D plot for the kurtosis of baryon number fluctuation κσ 2 as a function of the temperature and baryon chemical potential in the PNJL model. In Fig. 5 we show the phase transition lines and 2D plot for κσ 2 as a function of the baryon chemical potential along different freezeout lines for PNJL-1,PNJL-2 and PNJL-3, respectively. Different from the NJL model, in the PNJL model, when the gluodynamics is taken into account, there will be two phase transitions: one for the chiral phase transition and another for the deconfinement phase transition, and the deconfinement phase transition line lays below the chiral phase transition line at small chemical potential region. The two separate phase transition lines can be obviously seen from the (T, μ) phase diagram in Fig. 5. In the 3D plot Fig. 4 one can observe two separate phase boundaries at small baryon chemical potentials, and a valley forms in between the two phase boundaries. It is noticed that the magnitude of baryon number fluctuation κσ 2 is quite small along the chiral phase boundary, but around 1.5 along the deconfinement phase boundary.
By changing the coupling constant in the vector channel and parameter sets, the CEP of the chiral phase tran-sition in the PNJL model can shift. The CEPs are located at (μ E B = 919.1 MeV, T E = 123.5 MeV) and (μ E B = 979.5 MeV, T E = 104.2 MeV) in the PNJL-1 and PNJL-2 models, respectively. Even though the two critical baryon chemical potentials μ c B in the PNJL-1 and PNJL-2 are almost the same, the critical temperature in the PNJL-1 model is higher than that in the PNJL-2. For the parameters used in the PNJL-3, there is no CEP shows up in the phase diagram. From the 3D plot in Fig. 4, one can observe two obvious phase boundaries for the chiral and deconfinement phase transitions in the crossover side, and the magnitude of the deconfinement ridge is much higher than the chiral ridge. If a CEP for the chiral phase transition exists in the PNJL model, the structure of the CEP mountain for the chiral phase transition looks as the same as that in the NJL model, and one can observe a negative region of κσ 2 around the CEP above the chiral phase boundary extended from the crossover side. The only difference is that the deconfinement phase boundary extends to the CEP mountain and merges with the chiral phase boundary.
The structure of κσ 2 along the freeze-out line in the PNJL model is more complicated due to the two separated phase boundaries. Comparing with the NJL model, the magnitude of κσ 2 at small baryon chemical potentials in the PNJL model is in agreement with experiment measurement due to the contribution from gluodynamics. We also choose three different freeze-out lines defined in Eq. (18): (1) starting from the back ridge of the chiral phase boundary, goes through the negative region and then crosses the foot of the CEP mountain; (2) starting from the back ridge of the deconfinement phase boundary, and then crosses the foot of the CEP mountain; (3) starting from the deconfinement phase boundary and keeps far away from the CEP mountain. These three different freeze-out lines are indicated by long dashed, dashed-dotted and dashed lines, respectively. Here the back/ front ridge also means the higher/lower temperature side comparing with the phase boundary, respectively. The same as in the NJL model, the peak structure of κσ 2 along the freeze-out line in the PNJL model is solely related to the CEP mountain, when the freeze-out line crosses the CEP mountain or crosses the foot of the CEP mountain, there will be a peak showing up for κσ 2 , and the location of the peak is related to the location of the CEP mountain. If there is no CEP, κσ 2 keeps flat in almost the whole chemical potential region. In the PNJL-1 and PNJL-2 models with the existence of CEP in the phase diagram, it is found that for the first case of the freeze-out lines, one can observe the dip-peak structure for κσ 2 along the freeze-out lines. The freeze-out line starts from the back ridge of the chiral phase boundary and goes through the negative region of κσ 2 , then crosses the foot of the CEP mountain, therefore, κσ 2 decreases from around 0.5 and then down to negative value then rises up quickly and shows up a high peak around the critical chemical potential, thus shows a dip and then a peak structure of κσ 2 along the freeze-out line. For the second case, the freeze-out line starts from the back ridge of the deconfinement phase boundary, and has no chance to go through the negative region, then crosses the foot of the CEP mountain. For the third case when the freeze-out line is far away from the phase boundary, κσ 2 along the freeze-out line is almost flat.
In the PNJL-3 model, there is no CEP in the phase diagram, and no special structure of κσ 2 along the freeze-out lines is observed.

The kurtosis of the baryon number fluctuation κσ 2 in
the μPNJL model In Fig. 6, we show the 3D plot for the kurtosis of baryon number fluctuation κσ 2 as a function of the temperature and baryon chemical potential in the μPNJL model. In Fig. 7 we show the phase transition lines and 2D plot for κσ 2 as a function of the baryon chemical potential along different freeze-out lines for μPNJL-1 μPNJL-2 and μPNJL-3 models, respectively. Same as that in the PNJL model, in the μPNJL model, there also exist two separate phase transitions for the chiral phase transition and deconfinement phase transition, and the deconfinement phase transition line also lays below the chiral phase transition line at small chemical potential region. From the 3D plot Fig. 6 one can observe two separate phase boundaries at small baryon chemical potentials. The height of the ridge along the deconfinement phase boundary is around 1.5 at small chemical potentials in the μPNJL model, which is similar with in the PNJL model. Similar to that in the PNJL model, we choose three different freeze-out lines defined in Eq. (18): (1) starting from the back ridge of the chiral phase boundary, goes through the negative region and then crosses the foot of the CEP mountain; (2) starting from the back ridge of the deconfinement phase boundary, and then crosses the foot of the CEP moun-

The kurtosis of the baryon number fluctuation κσ 2 in a realistic PNJL model
In last section, we have investigated gluodynamics contribution to the baryon number fluctuations, and analyzed the formation of the dip and peak structures of the kurtosis along the imagined freeze-out lines. In this section, we will investigate the kurtosis along the experimental freeze-out line in a realistic 3-flavor PNJL model which takes into account 8quark interaction [45]. The effective potential is given below: where σ f = ψ f ψ f corresponds to quark condensates and f takes u, d for two light flavors while s for strange quark.
f with M f the dynamically generated constituent quark mass: If σ f = σ u , then σ f +1 = σ d and σ f +2 = σ s , and so on in a clockwise manner. U describes the contribution from self interaction of and¯ and it reads [46]: where and correspond to the effective potential of the Polyakov loop and the Jacobian of the transformation from the Polyakov loop to its trace, respectively. Besides, κ is a dimensionless parameter. b 2 (T ) is a temperature dependent coefficient which is chosen to have the form of   [42][43][44] are shown in dots, and the freeze-out temperatures and baryon number chemical potentials for lower energy heavy-ion collisions summarized in [47] are shown in squares, and two fitted freeze-out lines f 1, f 2 are shown by long dashed and dashed-dotted lines, respectively Follow [45], the parameters of the NJL part are fixed by vacuum properties and the parameters of Polyakov loop part are fixed by global fitting of the pressure density at zero chemical potential, and the details are listed in Tables 4 and  5, respectively. Fig. 10 The 2D plot for κσ 2 as a function of the baryon chemical potential (left) and the collision energy (right) in the realistic PNJL model along freeze-out lines f 1 (dashed line) and f 2 (dashed-dotted line) comparing with STAR Net-proton measurement [28]. The long dashed line is κσ 2 in a realistic NJL model [10]  With these parameters, as shown in [45], at zero chemical potential μ B = 0, the equation of state, baryon number fluctuations above the critical temperature are in good agreement with Lattice data. The kurtosis of baryon number fluctuation κσ 2 as a function of the temperature at zero baryon number density is shown in Fig. 8. It is noticed that κσ 2 in the realistic PNJL model in general is in good agreement with lattice data, especially comparing with the NJL model, PNJL model as well as μPNJL model. The phase transition line is shown in Fig. 9, and the CEP is located at (μ E B = 720 MeV, T E = 93 MeV). In Fig. 9, the freeze-out temperatures and baryon number chemical potentials extracted from BES-I at RHIC [44] are shown in dots, and the freeze-out temperatures and baryon number chemical potentials for lower energy heavyion collisions summarized in [47] are shown in squares, and the two fitted freeze-out lines f 1, f 2 described by T (μ) = 0.158 − 0.14μ 2 − 0.04μ 4 − 0.01 exp(−(μ − 0.067)/0.05) and T (μ) = 0.158 − 0.14μ 2 − 0.04μ 4 as used in [28] are shown in long dashed and dashed-dotted lines, respectively. Note that in these two formulas the T and μ B are in GeVs. The first freeze-out line f 1 starts from the back ridge of the phase boundary and f 2 starts from the front ridge of the phase boundary 9.
The kurtosis κσ 2 from the realistic PNJL model along the two freeze-out lines f 1, f 2 fitted from experimental data are shown in Fig. 10, the left figure is shown as a function of the baryon number chemical potential and the right figure is shown as a function of the collision energy, where we have used the following relation between the chemical potential and the collision energy: Note that in this formula the T and μ B are also in GeVs.
We can see that the kurtosis κσ 2 from the realistic PNJL model along the freeze-out line f 1, which starts from the back ridge of the phase boundary, develops a dip structure around μ B = 0.2 GeV( √ s = 20 GeV), and a peak structure at around μ B = 0.45 GeV ( √ s = 6 GeV). To our surprise, the kurtosis κσ 2 from the realistic PNJL model along this experimental freeze-out line agree with BES-I result very well. Along the second freeze-out line f 2, which starts from the front ridge of the phase boundary, the kurtosis only develops a peak structure at around μ B = 0.45 GeV ( √ s = 6 GeV) and no dip structure is developed. This supports our qualitative analysis in Sect. 2.

Conclusion and discussion
In this work, firstly we qualitatively investigate the kurtosis κσ 2 of net baryon number fluctuation and analyze the formation of its dip and peak structures along the imagined freezeout lines in the NJL model, PNJL model as well as μPNJL model with different parameter sets, and then we apply a realistic PNJL model and quantitatively investigate its κσ 2 along the real freeze-out line extracted from experiments.
Through qualitative analysis, we find that: (1) at zero chemical potential, the magnitude of κσ 2 is rather small in the NJL model comparing with lattice result, and it can reach around 1.5 in the PNJL and μPNJL models around the critical temperature. This indicates that gluodynamics plays important role in the baryon number fluctuation κσ 2 ; (2) the peak structure of κσ 2 along the freeze-out line is solely determined by the existence of the CEP mountain. When the freeze-out line crosses the CEP mountain or crosses the foot of the CEP mountain, there will be a peak showing up for κσ 2 along the freeze-out line, and the location of the peak is close to the location of the CEP mountain. The higher the peak, the closer the peak location to the CEP. In the case of no CEP, it is found that κσ 2 keeps flat almost in the whole chemical potential region and does not show any structure; (3) the formation of the dip structure is more complicated: firstly, it requires the existence of the CEP in the QCD phase diagram; secondly, it is sensitive to the relation between the freeze-out line and the phase boundary. The dip structure can be formed if the freeze-out line starts from the temperature above the critical temperature at μ = 0 and crosses the phase boundary from above when the CEP exists in the phase diagram. Because the CEP is for chiral phase transition and there is a negative region of κσ 2 from the crossover side, if the freeze-out line starts from the back ridge of the chiral phase boundary and if it goes through the negative region of κσ 2 , and then crosses the CEP mountain or the foot of the CEP mountain, in this case one can observe a dip structure and a peak structure, and the magnitude of κσ 2 will go to negative at the dip. If the freeze-out line starts from the back ridge of the deconfinement phase boundary, and then crosses the foot of the CEP mountain, we can also see a dip and peak structure of κσ 2 along the freeze-out line, but the magnitude of κσ 2 will not go to negative at the dip. Therefore we can read the information on how does the freeze-out line crosses the chiral phase boundary from the negative/positive value at the dip of measured κσ 2 along the freeze-out line.
Quantitatively, we use a reparameterized realistic PNJL model, with its critical temperature, equation of state and baryon number fluctuations in good agreement lattice data at zero chemical potential. To our surprise, the kurtosis κσ 2 produced from the realistic PNJL model along the experimental freeze-out line agrees with BES-I data well. This may indicate that the equilibrium result can explain the BES-I data on baryon number fluctuations. Indeed, from the analysis in [48], after collision, the system reaches thermalization quickly in quite high temperature and then evolves in equilibrium state, e.g. in the collision energy of √ s = 200 GeV, the system reaches thermalization at around T 210−230 MeV , which is much higher than the freeze-out temperature as well as the phase transition temperature. It is worth to point out that the extracted freeze-out temperatures from beam energy scan measurement are indeed higher than the critical temperatures at small chemical potentials, which supports our qualitative analysis on the formation of dip structure of κσ 2 along the freeze-out line.
At last, we should mention that in this work, even though our quantitative result from static thermodynamics can describe BES-I data [28] well, the non-thermal effect, or memory effect [49][50][51], and finite size effect deserves further studies to locate the CEP.