QCD chiral phase transition and critical exponents within the nonextensive Polyakov-Nambu-Jona-Lasinio model

We present a nonextensive version of the Polyakov-Nambu-Jona-Lasinio model that is based on nonextentive statistical mechanics. This new statistics model is characterized by a dimensionless nonextensivity parameter q that accounts for all possible effects violating the assumptions of the Boltzmann-Gibbs (BG) statistics (for , it returns to the BG case). Based on the nonextensive Polyakov-Nambu-Jona-Lasinio model, we discussed the influence of nonextensive effects on the curvature of the phase diagram at and especially on the location of the critical end point (CEP). A new and interesting phenomenon we found is that with an increase in q, the CEP position initially shifts toward the direction of larger chemical potential and lower temperature. However, when q is larger than a critical value , the CEP position moves in the opposite direction. In other words, as q increases, the CEP position moves in the direction of smaller chemical potential and higher temperature. This U-turn phenomenon may be important for the search of CEP in relativistic heavy-ion collisions, in which the validity of BG statistics is questionable due to strong fluctuations and long-range correlations, and nonextensive effects begin to manifest themselves. In addition, we calculated the influence of the nonextensive effects on the critical exponents and found that they remain almost constant with q.


I. INTRODUCTION
The QCD phase diagram is extremely important for us to be able to understand the evolution of strongly interacting matter, especially for our understanding of compact stars and the early universe [1][2][3][4]. When studying a QCD phase diagram, a statistical method based on Boltzmann-Gibbs (BG) statistics is often used. However, strictly speaking, BG statistics can only be applied to systems in equilibrium and within the thermodynamic limit. Obviously, in relativistic heavy-ion collisions, quarkgluon plasma (QGP) is produced, exhibiting strong intrinsic fluctuations and long-range correlations. In addition, the volume of QGP is small and evolves rapidly. Therefore, this system is far from uniform, and it is impossible to establish a global equilibrium [5]. As a result, the application of the common BG statistics in such collisions is questionable.
Thus, nonextensive statistics were first proposed by Tsallis, also known as Tsallis statistics [6]. The most typical feature of Tsallis statistics is the replacement of the usual exponential distribution factors by their q-exponential equivalents [7][8][9], where and its inverse function is (3) All factors that may violate BG statistical assumptions are summarized in the nonextensivity parameter q. For , Tsallis statistics returns to BG statistics because of and . Tsallis distribution has been observed in nature in several experimental and model systems. Impressive experimental examples include the following: (i) in high-energy physics, use of the Tsallis distribution to describe the transverse momentum spectra is now a standard practice [10][11][12][13][14][15][16][17][18][19][20], (ii) cold atoms in dissipative optical lattices exhibit an unusual transport behaviour that cannot be described within BG statistics [21], (iii) in Ref. [22] a particular case of the nonextensive scaling law in confined granular media is experimentally validated. From among the model systems, Refs. [23,24] show in a very clear manner that due to a break in ergodicity, the system has crossed from BG statistics to Tsallis statistics. In addition, a study of the nonextensive behavior of the QCD strong coupling constant is described in Ref. [25], in which the inconsistency between theory and experiment in the non-perturbative region is successfully considered. Finally, more about Tsallis statistics and its diverse applications can be found in Ref. [26].
Therefore, considering a real experimental environment in relativistic heavy-ion collisions, in this paper, we use Tsallis statistics to study the QCD phase transition. For this purpose, we generalize the PNJL model to its nonextensive version. Compared with the NJL model, this model can simulate the quark confinement effect and fit lattice QCD data more successfully by introducing a Polyakov-loop potential [27]. Moreover, other models, such as the linear sigma model and NJL model have also been generalized to the nonextensive version for studying the QCD phase diagram [8,9]. This paper is organized as follows: In Sec. II, we introduce the nonextensive version of the PNJL model. In Sec. III, we discuss the influence of nonextensive effects on the QCD phase transition, especially on the position of CEP and the critical exponents. Finally, we give a brief summary of our work in Sec. IV.

A. PNJL model
As a first step, we give a brief introduction to the PN-JL model. The Lagrangian of the two-flavor and threecolor PNJL model is as follows [27]: where represents the two flavor quark field with three colors, and with τ a (a = 1, 2, 3) stands for the current quark mass matrix. corresponds to the Pauli matrices in flavor space, and G is the effective coupling strength of the four point interaction of quark fields.
The effective Polyakov-loop potential is expressed in terms of the traced Polyakov-Loop expectation value and its conjugate (5) Index c refers to the color of the quark and . The Polyakov-loop L is defined as where is the temporal component of the Euclidian gauge field , , and P denotes the path ordering. Moreover, for simplicity, we take the approximation following Refs. [28][29][30][31][32][33]. According to Eq. (5), we have . The thermodynamic potential density function can be determined in the mean field approximation as where M is the effective quark mass and relates to the quark chiral condensate as and in which is the single quasi-particle energy, and is the quark chemical potential. In the above integrals, the cut-off only acts on the vacuum integral, whereas the medium dependent integrals have been extended to infinity [27,28,34,35]. U For the Polyakov-loop effective potential , we take the following two commonly used forms, for which the determination of the parameters is used to fit the pure gauge lattice data [27,29].
(1) The polynomial effective Polyakov-loop potential is [27,36,37] with a temperature-dependent coefficient and the corresponding parameters are given in Table 1.
(2) The Logarithmic effective Polyakov-loop potential is [29] with the temperature-dependent coefficients Φ,Φ ⩽ 1 the corresponding parameters are given in Table 2. Here, the logarithmic form constrains .
In a pure gauge sector, . However, the value of is adjusted to account for the presence of dynamical quarks. Here, we let following Ref. [38]. Moreover, the parameters for the NJL model part of the effective Lagrangian are summarized in Table 3, which are determined by fitting the pion decay constant and the pion mass 139.3 MeV [27].

Ω Φ
Finally, the solutions of the mean field equations are obtained by minimizing the thermodynamic potential function with respect to M and , that is B. Nonextensive PNJL model In short, the introduction of nonextensivity to a given model means that we need to make substitutions, as shown in Eq. (1) [8,9,39]. Moreover, we will use two simplifications in the following calculations, as in Ref. [40].
(i) The nonextensive effects are not considered in the pure Yang-Mills sector. That is, the Polyakov-loop potential remains unchanged and includes nonextensive effects implicitly only through the saddle point equations. µ (ii) The usual PNJL model parameters remain unchanged. We treat q just as researchers treat volume V in the study of finite-size effects, as a thermodynamic variable on the same footing as T and [41,42]. In fact, this is all based on the ansatz that the parameters determined at zero temperature and zero quark chemical potential can be used to study the finite temperature and finite quark chemical potential.
Therefore, within the nonextensive PNJL model, the thermodynamic potential density function becomes where Because the typical value of the nonextensivity parameter q is found to be in high-energy collisions [12,13,43,44], we only consider the case of in this paper.
In order to ensure that is always a non-negative real function, the following condition must be supplemented: If this condition is not met, for , one can use an approach with the Tsallis cut-off prescription For , the Tsallis cut-off prescription is where for particles and for antiparticles. Else, taking as an example, one can use an approach without the Tsallis cut-off prescription As pointed out in Ref. [8], there are certain problems with or without the Tsallis cut-off prescription. The Tsallis cut-off prescription limits the allowed phase space considerably. Without the Tsallis cut-off prescription the entropy makes a jump on the Fermi surface, which has not been observed in nuclear matter. Therefore, in order to avoid disputes about which cut-off should be used and to ensure the reliability of the results, Eq. (18) is always satisfied in our numerical calculations. Moreover, it is important to realize that for one always gets , as long as . This means that nonextensive effects only appear when the temperature is high enough.
In order to study the QCD phase diagram, based on Eq. (15), we need to solve the following coupled equations: where the q-version of the Fermi-Dirac distribution is and q → 1 For , they return to the distribution function of the usual PNJL model.

A. QCD phase transition
We plot M and as a function of T for four different q, ( , 1.05, 1.1, 1.15) as well as two different ( , ), as shown in Figs. 1 and 2. We found that even considering the nonextensive effect, the finite-temperature QCD transition is still not a real phase transition but a crossover [45]. However, as q increases, the transition occurs at a smaller pseudo-critical temperature . The same conclusion also appears in the nonextensive linear sigma model [9]. In addition, it is worth mentioning that the nonextensive effect has no influence on the QCD phase transition at zero temperature and finite chemical potential, as mentioned above.
κ For small values of the chemical potential, the shape of the QCD crossover transition line can be characterized by its curvature, which is accessible to lattice QCD [46,47]. Therefore, next we discuss the influence of nonextensive effects on the curvature . This is defined according to the Taylor series [48,49]: and herein, we choose . The pseudo-critical temperature is determined by the peak of the thermal susceptibility, as follows [33]: For the Polyakov-loop potential , the change in the crossover transition line with parameter q is shown in Fig. 3. We find that the pseudo-critical line gets lower and flatter with increasing q. In fact, the conclusion is the same for the Polyakov-loop potential . In order to better illustrate this, the change in curvature with q is listed in Table 4. Apparently, decreases as q increases. However, it should be noted that the curvature values we calculated are significantly higher than the values obtained in lattice QCD studies [46,47].  From  Fig. 4, we clearly see that when , the chiral condensation is discontinuous, which corresponds to a first-order phase transition. Moreover, when , is continuous, which corresponds to a crossover transition. Therefore, there must be a critical termination point for the first-order phase transition line in the middle chemical potential region, which is CEP. As we know, in the neighborhood of the CEP position, the susceptibility tends to diverge. Taking and the Polyakov-loop potential as an example, from Fig. 5,   we can clearly see that when , the susceptibility shows quite a sharp and narrow divergent peak, which determines the position of CEP as .
The influence of the nonextensive effect on the position of CEP is shown in Figs. 6 and 7. We found that its impact on the position of CEP is not as simple as previously thought. As q increases, the CEP position moves in the direction of the larger chemical potential and lower temperature at first, but then, when q is larger than a critical value , the CEP position moves in the opposite direction. In other words, as q increases, CEP moves in the direction of the smaller chemical potential and higher temperature. Obviously, this interesting U-turn phenomenon is independent of the choice of Polyakov-loop potentials. For and , the critical values are 1.1 and 1.08, respectively. It should be noted that this phenomenon was not found in Refs. [8,9]. Excluding the use of different models, one possible reason is that in Refs. [8,9], q was calculated only up to 1.1, which is just at the turning point. Our results are meaningful for the search of the CEP position in relativistic heavy-ion collision experiments. Based on our results, the search for areas with larger or smaller chemical potential depends on the critical value .

B. Critical exponents
As we all know, in the vicinity of CEP, the divergence in susceptibility can be described by the critical exponents. Regarding the critical exponents, there are two important physical concepts: the scale hypothesis and the universality assumption. Regarding the scale hypothesis, its basic idea is that when approaching the critical point, the correlation length , and the singularity of determines the singularity of all thermodynamic functions. From this, the scaling law that should be satisfied between the critical exponents can be derived. Regarding the universality assumption, it refers to a system with the same spatial dimension d and order parameter dimension n, with the same critical exponent and belonging to the same universal category. However, it should be pointed out that the validity of these concepts is based on the equilibrium phase transition system described by BG statistics. Therefore, in a system that deviates from the description of BG statistics, the critical exponents may not be completely determined by d and n, and the scaling law may need to be reconstructed or modified [50][51][52][53][54]. Based on this, in this subsection, we use Tsallis statistics to study the critical exponents and discuss the influence of the nonextensivity parameter q on them. However, it should be pointed out that the critical exponents calculated below are mean-field values, due to the mean-field approximation employed in this work.
Here, we choose a specific direction, which is denoted by ( ), for calculating the critical exponents using the path from lower (higher) T toward (represents ) with the quark chemical potential fixed at (represents ). Using the linear logarithmic fit, we obtain is the critical exponent of susceptibility, while is the critical exponent of the order parameter O, and , are constants. At and in the direction , the fitting procedure of the critical exponents for thermal susceptibility and quark mass is shown in Figs. 8 and 9.
The variation in critical exponents with q is shown in Tables 5 and 6. We found that when q increases from 1.0 to 1.1, the critical exponent remains almost unchanged, regardless of the Polyakov-loop potentials selected. However, for the critical exponent , taking the Polyakov-loop potential as an example, we found that for the direction ( ), decreases (increases) with an increase in q. For the Polyakov-loop potential , thisβ trend is the opposite. However, if we take the average value as the critical exponent parallel to the T axis, we find that is stable around 0.337 and hardly changes with q. In other words, even in Tsallis statistics, the critical exponents and are still the same as their mean-field values and in BG statistics [55][56][57][58].
It is worth noting that, in Ref. [50], the critical behavior of a two-dimensional Ising model with nonextensive statistics was studied, and it was found that for , the critical exponents , , and depend on q in the form , , and . However the critical exponent does not depend on q. Using the same model, in Ref. [51], it was found that although the critical expo-   nent is different from its value in BG statistics, it does not depend on q within the margin of error. In particular, the critical exponent depends on q in a linear manner. At present, the impact of the nonextensive effect on the critical exponent, especially for the two cases of and , is not very clear. Therefore, it is worthy of further research.

IV. SUMMARY AND CONCLUSION
In this paper, combined with the Tsallis statistics and the PNJL model, we investigated the sensitivity of the QCD phase transition and critical exponents with regard to deviations from usual BG statistics. Regarding the QCD phase diagram, we found that the influence of nonextensive effects on the CEP position shows a very interesting U-turn phenomenon. At the beginning, with an increase in q, the CEP position moves toward the direction of larger chemical potential and lower temperature. However, when q is larger than a critical value , as q increases, the CEP position moves in the opposite direction, that is, in the direction of smaller chemical potential and higher temperature. Because of this U-turn phenomenon, we found that searching for CEP in larger or smaller chemical potential regions in relativistic heavyion collisions depends on a critical value . This is a very interesting result that deserves further study. Regarding the critical exponents, numerical results based on Tsallis statistics show that the critical exponents remain almost constant with q. In other words, different from the findings in Refs. [50,51], we found that the critical exponents do not depend on the choice of Tsallis statistics or BG statistics. In addition, quark stars, as candidates for observed massive stars ( ), have attracted much attention in astronomy [4,[59][60][61]. Therefore, studying the influence of nonextensive effects on the structure and evolution of protoquark stars is a very meaningful topic. We shall study these issues in future.