Baryon number fluctuations and the phase structure in the PNJL model

We investigate the kurtosis and skewness of net-baryon number fluctuations in the Polyakov loop extended Nambu–Jona-Lasinio (PNJL) model, and discuss the relations between fluctuation distributions and the phase structure of quark-gluon matter. The calculation shows that the traces of chiral and deconfinement transitions can be effectively reflected by the kurtosis and skewness of net-baryon number fluctuations not only in the critical region but also in the crossover region. The contour plot of baryon number kurtosis derived in the PNJL model can qualitatively explain the behavior of net-proton number kurtosis in the STAR beam energy scan experiments. Moreover, the three-dimensional presentations of the kurtosis and skewness in this study are helpful to understand the relations between baryon number fluctuations and QCD phase structure.


Introduction
The exploration of QCD phase diagram and search for phase transition signatures of strongly interacting matter are subjects of great interest in high energy nuclear physics. Intensive searches on relativistic heavy-ion collision (HIC) have been performed in laboratories. The experimental data indicate that the transformation from quark-gluon plasma (QGP) to hadrons at small chemical potential and high temperature is a smooth crossover [1], consistent with the calculations of lattice QCD (LQCD) [2][3][4][5][6][7]. However, the QCD phase structure at large chemical potential is still unavailable in LQCD due to the sign problem in calculation.
The most fundamental issues in the non-perturbed region of QCD are the chiral symmetry breaking and quark confinement. Some QCD based models (e.g., [8][9][10][11][12][13][14][15][16]) and the Dyson-Schwinger equation approach [17][18][19][20][21] show that QGP undergoes a first-order chiral transition at large chemical potential. In the endpoint of the first-order transition there exists a crita e-mail: gyshao@mail.xjtu.edu.cn ical point (CP) connecting with the chiral crossover separation line. The important questions people concern are that whether the phase transitions can leave some traces at the final state in HIC experiments, and what the signatures are in measurements.
In theory, the particles fluctuations are most powerful and unique to investigate the properties of a thermal medium. Therefore, the event-by-event measurements of statistical distributions of fluctuations of conserved charges are crucial in searching for QCD phase structure [22,23]. The measurement of high-order fluctuations have been already suggested to explore the QCD phase transition (e.g., [24][25][26][27][28][29][30][31][32][33][34][35]). Recently, the preliminary result of the first phase of STAR Beam Energy Scan (BES-I) in Au+Au collisions indicates that the kurtosis of net-proton number fluctuation presents a non-monotonic energy dependence and a large deviation from the Poisson baseline [36][37][38]. This phenomenon can be qualitatively explained in the σ field model with a critical fluctuation of a first-order phase transition [27]. It possibly means that the experiments in Au+Au collisions with √ s = 7.7-27 GeV pass by the critical region of the firstorder chiral transition.
The measurements of net-proton fluctuation distributions and fluctuations of electric charge and strangeness greatly promote the investigation on the critical phenomenon of QCD phase transition. The baryon number fluctuations are investigated recently in different models, for example, the NJL model [39,40], the Polyakov Loop extended Quark Meson model [41], the functional renormalization group approach [42], the DSE approach [18], the interacting hadronic model [43,44] and the holographic QCD model [45]. Up to now the interpretation of the experimental results remains controversial [28][29][30][31]39,41,43,44] and phase transformation mechanism leading to the non-monotonic energy dependence of net proton kurtosis is still not clear.
In this study we will investigate the baryon number fluctuation distributions up to fourth order and discuss how they are affected by the chiral and deconfinement transitions in the PNJL model. The calculation shows that both the chiral and deconfinement transitions are related to the statistical distributions of net-baryon number fluctuation not only in the critical region but also in the crossover region. It means that the QCD phase structure can be deduced to a certain degree from the kurtosis and skewness of net baryon number fluctuation distributions. We also for the first time display the three-dimensional diagrams of kurtosis and skewness of net-baryon fluctuations in the PNJL model. These diagrammatic presentations are very helpful to understand the relations between baryon number fluctuation distributions and QCD phase structure.
The paper is organized as follows. In Sect. 2, we introduce the thermodynamical description of fluctuations of conserved charges in a thermal system and the 2 + 1 flavor PNJL quark model. In Sect. 3, we present the numerical results of the QCD phase diagram and the statistical distributions of net baryon number fluctuations. We then discuss the relation between fluctuation distributions and QCD phase structure, as well as the experimental results of net-proton number distributions. A summary is finally given in Sect. 4.

Fluctuations of conserved charges
For a grand-canonical ensemble in thermal equilibrium, the pressure of the system is related to the logarithm of the partition function [46]: where V and T are the volume and temperature of the system. The μ B , μ Q , μ S are the chemical potentials of conserved charges, the baryon number, electric charge and strangeness, respectively. The generalized susceptibilities of conserved charges can be derived by taking the partial derivatives of the pressure with respect to the corresponding chemical potentials [38] χ The cumulants of multiplicity distributions of the conserved charges are connected with the generalized susceptibilities by In experiments, observables are constructed by the ratio of cumulants, which cancel the volume dependence and then can be compared with theoretical calculations of the generalized susceptibilities.
For an arbitrary distribution we can define the Gaussian width σ 2 with the non-Gaussian fluctuations. In this study, we focus on two important statistic quantities, the skewness (Sσ ) and Kurtosis (κσ 2 ). For the net baryon number fluctuations, the skewness and Kurtosis are related to the high order cumulants as Similar relations can be derived for the fluctuations of electric charge and strangeness. We will calculate the susceptibilities of conserved charges in the three-flavor PNJL model. The chemical potentials μ B , μ Q , μ S used in experiments and LQCD simulations are related to the quark chemical potentials with the following relations and where μ u,d,s are the quark chemical potentials for up, down and strange quarks. The Lagrangian density in the 2 + 1 flavor PNJL model is taken as where q denotes the quark fields with three flavors, u, d, and s;m 0 = diag(m u , m d , m s ) in flavor space; G and K are the four-point and six-point interacting constants, respectively. Theμ = diag(μ u , μ d , μ s ) are the quark chemical potentials which are related to chemical potentials of the conserved charges through Eqs. (5) and (6). The covariant derivative in the Lagrangian is defined as where β = 1/T is the inverse of temperature and A 4 = i A 0 . The Polyakov-loop effective potential given in [47] is where The parameters a i , b i listed in Table 1 are fitted according to the lattice simulation of QCD thermodynamics in pure gauge sector. And T 0 is found to be 270 MeV as the critical temperature for the deconfinement phase transition of gluon part at zero chemical potential [48]. When fermion fields are included, a rescaling of T 0 = 210 MeV is implemented.
In the mean field approximation, the constituent quark mass can be derived as where φ i stands for quark condensate of the flavor i.
The thermodynamical potential is derived as where All the thermodynamic quantities relevant to the bulk properties of quark matter can be obtained from Ω. In the calculation a cut-off Λ is implemented in 3-momentum space for divergent integrations. We take the model parameters obtained in [49]: Λ = 603.

Numerical results and discussions
In this section, we first analyze the QCD phase structure, and then investigate the statistical distributions of baryon number fluctuations and discuss their relations with the chiral and deconfinement transitions. In the calculation, we take μ Q = μ S = 0 for simplicity, which is approximately consistent with the experimental measurements [50]. Strange quark chemical potential is set to zero, which approximatively fulfills the requirement of strangeness conservation in the process of strong interaction.

Phase structure in the PNJL model
Since there have been already some researches on the phase structure in the PNJL model in literature (e.g. [8][9][10][11]), for the convenience of later discussion, here we just focus on some respects closely related to the baryon number fluctuations.
We show the quark condensate as a function of temperature for different quark chemical potentials in the upper panel of Fig. 1. It indicates that φ l /φ 0 decreases continuously when quark chemical potential μ q is equal to 0, 100, 200 and 250 MeV, respectively. This means the chiral phase transition is a smooth crossover at small chemical potential. For the case of μ q = 320 MeV, there is a jump of chiral condensate, indicating the appearance of the first-order transition.
The lower panel of Fig. 1 shows the partial derivatives of φ l respect to temperature. The evident features for μ q = 0 is that ∂φ l /∂ T has two maxima. This structure only exists for small chemical potential. The smaller maximum on the lower-temperature side will gradually disappear with the increase of quark chemical potential, as shown in the curves of μ q = 100 and 150 MeV, and finally only a single peak exists. The chiral crossover separation line can be derived by  connecting the peaks of the larger maxima in the T − μ q plane, which will be plotted later in the full QCD phase diagram. We also plot the three-dimensional diagram of ∂φ l /∂ T in Fig. 2. This figure clearly demonstrates that how fast the chiral crossover transition takes place as a function of temperature and chemical potential. It also indicates the existence of the first-order phase transition at large chemical potentials where ∂φ l /∂ T diverges.
We present the value of Φ and its derivative respect to temperature in Fig. 3. This figure shows that Φ and ∂Φ/∂T changes continuously for μ = 0, 100, 150, 200 and 250 MeV. There is a jump for μ = 320 MeV, which is induced by the restoration of chiral symmetry accompanied with the first-order chiral transition. The deconfinement transition line can be derived by the requirement that ∂Φ/∂T takes the maximum value.
However, there is a special case that we needs to pay attention to. Along the first-order transition line, although Φ and ∂Φ/∂T are discontinuous for μ q much larger than the critical chemical potential of chiral transition, the deconfinement still does not occur. The main reason is that the values of Φ are quite small on the both sides of the chiral first-order transition at low temperature. This forms the quarkyonic phase in which the chiral symmetry is restored but quarks are still confined. As a matter of fact, there also exists the other local maximum of ∂Φ/∂T at a higher temperature above the firstorder transition, as shown in Fig. 3 with μ = 320 MeV. If we connect these points with those obtained in the crossover region with ∂Φ/∂T taking the local maximum, we can obtain the deconfinement transition line in the full QCD phase diagram, as plotted in Fig. 5. The derivative of Φ respect to T in the lower panel of Fig. 3 reflects the width of the deconfinement transition. We visualize ∂Φ/∂T as a function of temperature and chemical potential in the three-dimensional diagram in Fig. 4, which clearly presents the features of the deconfinement transition. We notice that the peak of ∂Φ/∂T above the first-order transition line becomes quite flat. Figures 1, 2, 3 and 4 indicate that the chiral and deconfinement transitions at small chemical potentials occur in a wide range of temperature. To indicate the relation between chiral and deconfinement transition, we plot the complete phase diagram in Fig. 5. The black solid line is the chiral first-order transition line, and the black dash line is that of the chiral crossover separation line. The red dash curve is the deconfinement transition line obtained with ∂Φ/∂T taking The area filled with oblique lines is roughly the region of chiral crossover, and the yellow band is region of rapid chiral crossover (corresponding to the green peak in Fig. 2). For the deconfinement transition, the rapid transition region lies in the blue band. Figure 5 indicates that the chiral and deconfinement transitions can be approximately separated in the crossover region, but they almost overlap near the critical region. Comparing with the chiral transition, we find the blue band of the rapid deconfinement almost overlaps with the smaller peak of chiral transition, as shown in Fig. 2. This means the rapid deconfinement leads to the partial restoration of chiral symmetry at small chemical potential, but it is not strong enough to force the complete restoration of chiral symmetry. With the increase of chemical potential, it becomes difficult to separate the two kinds of phase transitions, in particular, in the critical region. For the spinodal structure (the gray area in Fig. 5) of chiral first-order transition, one can refer to Ref. [11] for a detailed description.

Kurtosis and skewness of baryon number fluctuations
In this subsection, we investigate the kurtosis and skewness of net-baryon number fluctuations and discuss their relations with the chiral and deconfinement transitions.
We show the contour map and three-dimensional landscape of the kurtosis κσ 2 of net-baryon number fluctuation in Figs. 6 and 7, respectively. In Fig. 6, the values of κσ 2 are negative in the red area and positive in the rest area. This figure has the same topological structure as the contour map derived in the σ field model [27] and the NJL model [39,40]. However, the main difference in the PNJL model is that there are two branches of the net-baryon number kurtosis in the contour map, and they converge in the critical region. One branch is along the chiral transition line, as derived in the NJL Fig. 6 Contour lines of the net-baryon kurtosis in the T − μ q plane. κσ 2 is negative in the red area and positive in the rest area Fig. 7 Three-dimensional structure of net-baryon kurtosis as a function of temperature and quark chemical potential model (Only this branch exists in the NJL model) [39,40]. For this branch, the value of κσ 2 is universally negative in a narrow region (red area in Fig. 6) when the critical point is approached on the crossover side of the chiral transition line. The higher order fluctuations become more and more intensive in the vicinity of the critical point. Figure 6 shows that the second branch of κσ 2 of netbaryon number fluctuation appears only in the crossover region. The peaks of the contour lines of this branch lie in the region of rapid deconfinement transition, as shown in Fig. 5. Such a structure can also be seen in the threedimensional diagram in Fig. 7. The appearance of this branch can be attributed to the partial restoration of chiral symmetry induced by the rapid deconfinement transition, as discussed in the last subsection.
The existence of two branches of κσ 2 indicates that both the chiral and deconfinement transitions are important to the density fluctuations. It is not difficult to understand such a structure in the PNJL model. Different from the NJL model and the σ field model, in the PNJL model there are two kinds of order parameters, quark condensates and Polyakov loop in the quark distribution function (closely related to the density fluctuations). The separation of the chiral and deconfinement transitions leads to the appearance of the branch in the crossover region. Figures 6 and 7 indicate that the statistical distributions of baryon number fluctuation are related to both the chiral and deconfinement transitions, which can be taken to explore the QCD phase structure in experiments.
The numerical results in Figs. 6 and 7 show that the value of κσ 2 deviates from unity in a wide region. κσ 2 increases quickly in the critical region and diverges at the critical point. It indicates that the magnitude of deviation from unity of κσ 2 depends on the distance from the critical point. If the chemical freeze-out curve passes by the critical region, even not very close to the critical point, a large deviation of baryon number kurtosis from unity possibly appears. The contour of net-baryon number kurtosis derived in the PNJL model can qualitatively explain the deviation from the Poisson baseline and the non-monotonic behavior of net-proton number kurtosis [36][37][38]. On the other hand, we note that the actual location of the QCD critical point is currently unknown due to the limitation of LQCD at large chemical potential. However, the topological structure of kurtosis of density fluctuation is universal for a system with a first-order phase transition and a critical point [51]. As for how the critical fluctuations propagates, it involves the dynamics of fluctuation in an expanding system [52][53][54]. More investigations on hydrodynamics and non-equilibrium effect are needed for the critical fluctuations [55,56].
Finally, we demonstrate the numerical results of the skewness Sσ of net baryon number fluctuation in Figs. 8 and 9. Figure 8 indicates that the value of Sσ is negative in a slender region close to the chiral transition line on the right side, and the Sσ increases quickly in the critical region on the left side of the chiral transition line. Figure 9 demonstrates the three-dimensional structure of the skewness as a function of temperature and chemical potential. From this figure it is easy to understand the increase of net-proton Sσ in experiments with the decrease of collision energy [36][37][38]. Similar to the distribution of κσ 2 , in the crossover region the partial restoration of chiral symmetry induced by the rapid deconfinement also affects the distribution of Sσ . In general, the value of Sσ is smaller than κσ 2 , therefore the measurement of κσ 2 in experiments is more effective to investigate the critical phenomenon of QCD phase transition.

Summary
We investigated the kurtosis and skewness of net-baryon number fluctuations within the PNJL model, and discussed their relations with the QCD phase structure. The calculation indicates that the QCD phase structure is mainly determined by the chiral and deconfinement transitions, as well as their coincidence and separation. It also shows that the statistical distributions of baryon number fluctuations are closely related the QCD phase structure.
Therefore, the measurement of the statistical distributions of particle fluctuations can provide important information about the QCD phase transition, in particular, the critical phenomenon. Compared with the preliminary result of the STAR beam energy scan, the calculation shows that the contour map of net-baryon number kurtosis in the PNJL model can qualitatively explain the deviation from the Poisson baseline of net-proton number kurtosis and the non-monotonic energy dependence. The three-dimensional presentations of baryon number multiplicity distributions in this study are helpful to understand the relations between baryon number fluctuation distributions and QCD phase structure. Further investigations with more fundamental theory are needed to give a quantitative description of the experimental data.