Thermodynamic properties and transport coefficients of QCD matter within the nonextensive Polyakov–Nambu–Jona-Lasinio model

We present a non-extensive version of the Polyakov-Nambu-Jona-Lasinio model which is based on the non-extentive statistical mechanics. This new statistics is characterized by a dimensionless non-extensivity parameter $q$ that accounts for all possible effects violating the assumptions of the Boltzmann-Gibbs statistics (when $q\rightarrow1$, it returns to the Boltzmann-Gibbs case). Using this q-Polyakov-Nambu-Jona-Lasinio model and including two different Polyakov-loop potentials, we discussed the influence of the parameter $q$ on chiral and deconfinement phase transition, various thermodynamic quantities and transport coefficients at finite temperature and zero quark chemical potential. We found that the Stefan-Boltzmann limit is actually related to the choice of statistics. For example, in the Tsallis statistics, the thermodynamic quantities $\frac{\epsilon}{T^{4}}$, $\frac{p}{T^{4}}$ and $\frac{s}{T^{3}}$ all increase with $q$, exceed their usual Stefan-Boltzmann limits and tend to a new $q$-related Tsallis limit at temperature high enough. Interestingly, however, due to a surprising cancellation, the high temperature limit of $c_{s}^{2}$ is still its SB limit $1/3$. In addition, we found some similarities between the non-extensive effect and the finite-size effect. For example, as $q$ increases (size decreases), the criticality of $\frac{c_{v}}{T^{3}}$ and $c_{s}^{2}$ gradually disappears. Besides, in order to better study the non-extensive effect, we defined a new susceptibility and calculated the response of thermodynamic quantities and transport coefficients to $q$. And found that their response patterns are different.


I. INTRODUCTION
Among all standard studies of the QCD matter, a statistical approach often used is Boltzmann-Gibbs (BG) statistics. However, strictly speaking, this approach is correct only when the corresponding heat bath is homogeneous and infinite. Obviously, in reality, this condition cannot always be met. Especially in the relativistic heavy-ion collisions, in which the quark-gluon plasma (QGP) produced experiences strong intrinsic fluctuations and long-range correlations. The size of QGP is small enough and it evolves rapidly. Therefore, this system is far from being uniform and no global equilibrium is established. As a result, some quantities become non-extensive and develop power-law tailed rather than exponential distributions. In such cases the application of the usual BG statistics is questionable.
Thus, a non-extensive statistics that extended BG statistics was first proposed by Tsallis [1]. The most typical feature of Tsallis statistics is that it replaces the usual exponential factors by their q-exponential equivalents [2][3][4], where exp q (x) = [1 + (1 − q)x] 1 1−q , (2) * zhaoyapeng2013@hotmail.com correspondingly, its inverse function is The non-extensivity parameter q represents all possible factors that do not satisfy the BG statistical assumptions. When q → 1, exp q (x) → exp(x), ln q (x) → ln(x) and Tsallis statistics returns to BG statistics. In high-energy physics, using Tsallis statistics to describe the transverse momentum distributions is now a standard practice [5][6][7][8][9][10]. It is excellent in meeting the experimental data, as pointed out by the PHENIX [11] and STAR [12] Collaborations at RHIC and by the ALICE [13], ATLAS [14] and CMS [15] Collaborations at the LHC. In addition, more and more physical branches, even biology, economics are described by Tsallis distribution. A general overview on Tsallis' statistics and its diverse applications can be found in Ref. [16]. Finally, it should be noted that Ref. [17] also raises doubts about the application of Tsallis statistics to relativistic heavy-ion physics.
Studying the thermodynamic properties and transport coefficients of the QCD matter has always been a matter of great interest to people. At high baryon density and low temperature, they are relevant to the study of compact stars [18][19][20][21]. For example, the equation of state can be combined with the Tolman-Oppenheimer-Volkoff (TOV) equation to study the mass-radius relationship, the internal structure of compact stars and further to study the tidal Love number k 2 and tidal deformability Λ. At high temperature and low baryon density, they are relevant to the QGP produced in relativistic heavy-ion collisions. Especially, a low value of the shear viscosity to entropy (η/s) is needed to explain the elliptic flow data [22], which means the QGP is actually a stronglycoupled medium. In addition, studies on thermodynamic quantities and transport coefficients may also help to reveal QCD phase transitions or a rapid crossover [23][24][25][26].
In this paper, the question we are concerned with is that when a strongly interacting system is described by Tsallis statistics, what is the difference between the thermodynamic quantities and transport coefficients and that of BG statistics. For this purpose, we generalize the PNJL model to its non-extensive version. Compared with NJL model, this model has proven to be more successful in reproducing lattice data concerning QCD thermodynamics [27]. Besides, other models such as the linear sigma model and NJL model have also been generalized to its non-extensive version to study the thermodynamic quantities of the QCD matter and its phase diagram [3,4]. This paper is organized as follows: In Sec. II, we introduce the non-extensive version of the PNJL model and discuss the q-dependence of the chiral and deconfinement phase transition at finite temperature and zero quark chemical potential. In Sec. III we analyze in detail the influence of the parameter q on the thermodynamic quantities and transport coefficients. Finally, we give a brief summary of our work in Sec. IV.

II. NON-EXTENSIVE PNJL MODEL:Q-PNJL
Before introducing the q-PNJL model, let's make a basic introduction to the PNJL model. The Lagrangian of the two-flavor and three-color PNJL model reads [27] where Ψ = (u, d) andm = diag(m u , m d ) with m µ = m d = m stands for the current quark mass matrix. τ a (a = 1, 2, 3) are Pauli matrices acting in flavor space and G is the effective coupling strength of four point interaction of quark fields. The effective Polyakov-loop potential U(Φ,Φ; T ) accounts for the self-interactions of the gauge field in which the normalized color-traced Polyakov-loop expectation value Φ and its Hermitian conjugationΦ are defined as where the Polyakov line is defined as and A 4 = iA 0 is the temporal component of Euclidian gauge field ( A, A 4 ), β = 1 T , and P denotes the path ordering. The covariant derivative is determined as here A µ = gA a µ λ a 2 and g is the SU (3) c gauge coupling. The λ a stands for the Gell-Mann matrices with λ 0 = 2 3 1. Under the mean-field approximation, the thermodynamic potential density function is where M means the dynamical quark mass. It relates to the quark chiral condensate σ = Ψ Ψ as follows and in which E p = p 2 + M 2 is the single quasi-particle energy. In the above integrals, following Refs. [27][28][29][30], the vacuum integral has a cut-off Λ whereas the medium dependent integrals have been extended to infinity.
Finally, the solutions of the mean field equations are obtained by minimizing the thermodynamic potential function Ω with respect to M , Φ andΦ, that is at vanishing chemical potential,Φ = Φ.

A. Polyakov-loop potentials
The functional form of the effective Polyakov-loop potential U that can be constructed from the center symmetry of the pure gauge sector is not unique. The required parameters are based on the pure gauge lattice data. Next two effective Polyakov-loop potentials are introduced.
The polynomial effective Polyakov-loop potential is [27,31,32]  with a temperature-dependent coefficient and the corresponding parameters are given in Table I.
In a pure gauge sector, T 0 = 270 MeV. However, in the presence of dynamical quarks, the critical temperature T 0 will have an N f dependence T 0 (N f ). For massless flavors, T 0 (2) = 208 MeV with an uncertainty of about 30 MeV. If we consider that quark has mass, the critical temperature will be lower. Here, we let T 0 (2) = 192 MeV follows Ref. [33]. Besides, it should be noted that with this Polynomial potential U P , the Polyakov-loop expectation value Φ (at µ = 0, Φ =Φ) will be greater than one at a temperature of a few hundred MeV and when T → ∞, Φ ≃ 1.11. The Logarithmic effective Polyakov-loop potential is [34] with the temperature-dependent coefficients and the corresponding parameters are given in Table II. Here, the logarithmic form constrains Φ,Φ ≤ 1.
Besides, the parameters for the NJL model part of the effective Lagrangian L P N JL are summarized in Table III.

B. q-PNJL model
In short, when we use Tsallis statistics instead of BG statistics to describe a system, it means we need to do the replacement as shown in Eq. (1). In this first case study we shall however take up the following two simplifications: (i) in the present treatment non-extensive effects are not considered in the pure Yang-Mills sector.
As a consequence, the Polyakov-loop potential remains unchanged and feels non-extensive effects implicitly only through the saddle point equations.
(ii) we shall not use any modifications to the usual PNJL model parameters due to non-extensive effects. We treat q as a thermodynamic variable in the same footing as T and µ. Similarly, in the study of finite-size effects, they also treat volume V as a thermodynamic variable in the same footing as T and µ. Fitting the parameters at T = 0, µ = 0 and V = ∞ and then studying the finitesize effects at finite temperature or/and quark chemical potential [35,36]. 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. Of course, it is also pointed out in the Refs. [37,38] that the parameter of the coupling constant G should depend on the order parameter Φ or Ψ Ψ and then implicitly on the temperature and the quark chemical potential. But here, we do not consider this situation.
Thus, within the q-PNJL model, the thermodynamic potential density function becomes In order to ensure that e q (x) is always a non-negative real function, the following constraint must be met.
And in this paper, as a first step, we consider only q > 1. This is because on the one hand the typical value of the non-extensivity parameter q for high energy collisions is found to be 1 ≤ q ≤ 1.2 [7,8,39,40]. On the other hand in the case of q > 1, q − 1 is a measure of intrinsic fluctuations of the temperature in the system considered [41,42], whereas q < 1, the interpretation of q − 1 is inconsistent [43,44]. So, in the case of zero quark chemical potential and finite temperature, the condition Eq. (19) is naturally satisfied. If not, we can use the following Tsallis cut-off prescription for q > 1 or without Tsallis cut-off prescription It should be pointed out that it is not clear so far under which circumstances Tsallis cut-off or without Tsallis cut-off should be used. A more detailed discussion can be found in Ref. [3].
Besides, it is important to realize that for T → 0 one always gets Ω q → Ω as long as q > 1. This means that we can expect any non-extensive signature only for high enough temperatures.
For studying the phase transition within the q-PNJL model at zero quark chemical potential and finite temperature, according to Eq. (11), the coupled nonlinear equations for the M and Φ can then be obtained as follows where the q-version of the Fermi-Dirac distribution is andn when q → 1, they return to the distribution function of the usual PNJL model.
152 142 a function of T for three different q, (q = 1, 1.05, 1.1) as well as two different U (U P , U L ). First of all, the finite-temperature QCD transition is not a real phase transition, but a crossover as Ref. [45] shows, and it is independent of the q and U. Secondly, we find that the influence of q on the chiral transition and deconfinement transition is consistent. That is to say, as q increases, the transition occurs early, even though there will be obvious quantitative differences for two different U. The same conclusion also appears in the non-extensive NJL model [3] and the non-extensive linear sigma model [4].  In order to better study the influence of parameter q on the crossover transition, we introduce the susceptibility, which is defined as The peak position corresponds to the pseudo-critical temperature T χ of the chiral transition and T d of the deconfinement transition, respectively. The results for two different U are presented in Tables IV, V.  From Tables IV, V we can clearly see that as q increases, both T χ and T d decrease. Specifically, as q increases from 1 to 1.1, the pseudo-critical temperatures T χ and T d decrease by approximately 6% − 8.2% and 9.4% − 11.3%, respectively. Besides, we find T χ > T d and the same result can be seen in Ref [46]. At last, we define a new susceptibility to describe the response of X to the parameter q. It can be seen from Fig. 3 that for M , the maximum value of |χ Mq | appears near the pseudo-critical temperature T χ . That is to say, at T χ , M has the largest response to q. However, it should be pointed out that in the case of the effective potential U L , the situation is slightly different. That is, |χ Mq | will have maximum values near the pseudo-critical temperatures T d and T χ , respectively. Regarding χ Φq , it shows a maximum value near the pseudo-critical temperature T d for both effective potentials. As for the response of other thermodynamic quantities and transport coefficients to q, we will discuss in detail in Sec. III.

III. QCD THERMODYNAMIC QUANTITIES AND TRANSPORT COEFFICIENTS
In this section, we mainly study the influence of parameter q on the thermodynamic quantities and transport coefficients within q-PNJL model. The reason why we are interested in thermodynamic quantities and transport coefficients is because they are not only sensitive to phase transition but also they can offer important information to other fields, like hydrodynamical models of the QGP, cosmological models of the early universe and models of massive objects in astrophysics as we emphasized above.

A. QCD thermodynamic quantities
All the thermodynamic information of a system is contained in the grand canonical potential which is given  by Ω in Eq. (17) evaluated at the mean-field extent. The pressure is with the vacuum normalization p q (0) = 0. The entropy density s q and energy density ǫ q are defined as follows In the SB limit, the QCD pressure for N 2 c − 1 massless gluons and N f massless quarks is given by [33] p SB T 4 = (N 2 c − 1)  where the first term is the gluonic contribution and the second term is the quark's contribution. For effective potentials U P , U L , at zero quark chemical potential and N c = 3, N f = 2, we have p SB /T 4 ≃ 4.06. Correspondingly, and s SB /T 3 = p SB /T 4 + ǫ SB /T 4 ≃ 16.23.
From Figs. 4, 5, 6 we can see that when q = 1, p/T 4 , ǫ/T 4 and s/T 3 all tend to their SB limit. However, as q increases, they increase rapidly until they exceed their corresponding SB limits. Taking ǫ/T 4 as an example, for the Polyakov-loop potential U L and the temperature is  fixed at 350 MeV. When q = 1, the value of ǫ/T 4 is 11.3, very close to the SB limit 12.18. But when q = 1.1, the value of ǫ/T 4 is 18.1, which is increased by 60%. For p/T 4 and s/T 3 , this value is 67% and 61%, respectively. This is the result we expect because in the q-PNJL model, we use Tsallis statistics instead of BG statistics, and q − 1 describes exactly the deviation from BG statistics. If q−1 is larger, the deviation from the SB limit is larger. That is to say, when a system is described by Tsallis statistics, then its high temperature limit is not the SB limit but the q-dependent Tsallis limit. This can be understood from the fact that the thermodynamic potential density function Ω q = Ω even when T → ∞. That is to say, for the high temperature limit, the non-extensive effect still exists. This is different from the finite-size effect. Consider a cube system with a length of L, as T increases, the thermal de Broglie wavelength decreases, and the effective size of the system becomes larger. Therefore, when T L → ∞, we can still consider it as an ideal gas system and the related physical quantities tend to their SB limits. This means that the finite-size effect does not change the SB limit [35,36,47]. Furthermore, from Figs. 7, 8, we find that the response patterns of ǫ/T 4 and s/T 3 for q are almost the same. They reach a maximum near the pseudo-critical temperature T c , then decrease and tend to be stable. About p/T 4 , from Eq. (29), it is known that χ pq = χ sq − χ ǫq . The result is shown in Fig. 9, in which its response increases with temperature and gradually stabilizes. For relativistic heavy-ion collisions, the speed of sound is an important quantity, and its square at constant entropy is defined by where c vq denotes the specific heat at constant volume, defined as From Fig. 10, like p/T 4 , ǫ/T 4 and s/T 3 , when q = 1, c v /T 3 tends to the usual SB limit (48.69 for U P and U L ) with increasing temperature. However, as q increases, the high temperature limit value also increases. In order to better show the effect of the non-extensivity parameter q on c v /T 3 , it could be interesting to normalise all lines with respect to their high temperature limit, as shown in Fig. 11. For U P , at q = 1, the normalized c v /T 3 starts to rise with increasing temperature T , then reaches the maximum near the pseudo-critical temperature T χ of chiral transition, and eventually tends to the SB limit. But for U L , there are two peaks. The first corresponds to the pseudo-critical temperature T d of deconfinement transition, and the second, although not obvious, corresponds to the pseudo-critical temperature T χ of chiral transition as Refs. [26,36] show. Regarding the dependence of the normalized c v /T 3 on q, the height and the pseudo-critical temperature T c of its peak decreases as q increases for two Polyakov-loop potentials U. Especially for U L , this peak flattens as q increases so that it is difficult to show the pseudo-critical temperature T c of the crossover transition. That is to say, as q increases, the critical behavior of c v /T 3 is smoothed out and this phenomenon also appears in the c 2 s , which we will discuss next.
The behavior of c 2 s is shown in Fig. 12. For two different U and q = 1, near the pseudo-critical temperature T c it has a dip and then approaches the ideal gas value of 1/3 at high enough temperatures.
The same conclusion appears in other models, such as the NJL model and the Polyakov-Quark-Meson (PQM) model [46,48]. Besides, as the q increases, we can see that the dip is gradually disappearing. This is similar to Ref. [49] where it calculated the speed of sound as a function of temperature for different q-values for a hadron resonance gas and found that when q is larger than 1.13, all criticality disappears. Interestingly, similar phenomena have appeared in the study of finitesize effects. Refs. [35,36,50] indicate that as the size decreases, the critical behavior of c 2 s also gradually or even completely disappears. It is reasonable to agree with each other because the finite-size effect is part of the nonextensive effect. In addition, it is worth noting that due to a surprising cancellation, the high temperature limit of c 2 s is still its SB limit 1/3, independent of U. Taking U L as an example, we find that at high temperature T = 350 MeV, when q increases to 1.1, s/T 3 and c v /T 3 increase by 61.4% and 59.4%, respectively. As the temperature further increases, the growth rates of s/T 3 and c v /T 3 tend to be the same. Therefore, the high temperature limit of c 2 s = s T 3 / cv T 3 is not affected by non-extensive effects.
Another quantity that is related to c 2 s is the interaction measure ∆ q (T ) which measures the deviation from the equation of state of an ideal gas ǫ = 3p due to interactions and/or finite quark masses, defined as and it is related to the c 2 The reason for doing this approximation comes from Refs [46,51], from which we can see that c 2 s and p/ǫ are in good agreement at both low temperature (< T c ) and high temperature (> 2.5T c ). At the intermediate temperature (T c ∼ 2.5T c ), c 2 s is slightly larger than p/ǫ. Therefore, it can be known from Eq. will cause the maximum value of ∆ q (T ). At the high temperature limit, c 2 s tends to 1/3 and ∆ q (T ) tends to zero. From Fig. 13 we can clearly see that it has a peak near the pseudo-critical temperature T c and then tends to zero. And the peak moves with q towards a lower temperature, which is similar to c v /T 3 and independent of U.

B. Transport coefficients
In this subsection we are mainly concerned with the influence of nonextensivity parameter q on the transport coefficients, such as shear viscosity η, electrical conductivity σ and bulk viscosity ζ. Based on linear σ model [52], NJL model [53][54][55][56][57][58], PQM model [59], PNJL model [60] and the Parton-Hadron-String Dynamics (PHSD) transport approach [61], we get a gross summary about the temperature dependence of these transport coefficients. η(T ) and σ(T ) decrease with temperature increase in the hadronic phase, while they increase with temperature increase in the QGP phase, and show a minimum at the transition temperature. For certain materials, like helium, nitrogen and water, this temperature dependence of η has been experimentally confirmed [62]. While ζ(T ) follows an opposite trend, which shows a maximum at the transition temperature [52,54,58,63].
The mathematical expressions of transport coefficients calculated from relaxation time approximation (RTA) in the kinetic theory approach and calculated from the one-loop diagram approximation in the quasi-particle Kubo approach are equivalent, as follows [50,64] It should be noted that n q andn q are not the usual Fermi-Dirac distribution but the q-version of the Fermi-Dirac distribution, defined by Eqs. (24), (25). In addition, all our discussions are based on a constant value of relaxation time τ = 1/Γ = 1 fm. The changes of η and σ with T and three parameters q are shown in Figs. 14, 15. First of all, we find that η and σ rise monotonically with temperature increase as Refs. [50,64] show, and they increase as q increases. More specifically, for η , its response to q increases monotonically with temperature, as shown in Fig. 16. For σ, its response to q appears an extreme value near the pseudo-critical temperature and then rises monotonically, as shown in Fig. 17. Regarding ζ, from Fig. 18 we find that the situation will be different for different U. For example, in the effective potential U P , ζ increases with temperature and has a significant maximum near the pseudo-critical temperature T c , and then tends to zero. However, in the effective potential U L , ζ will have two maximum values as the temperature rises. The first maximum value is not so obvious, corresponding to the vicinity of the pseudocritical temperature T d of deconfinement transition, and the second maximum value corresponds to the vicinity of the pseudo-critical temperature T χ of chiral transition. A similar double-peak structure also appears in the NJL model [53] and PNJL model [50]. Moreover, Ref. [50] indicates that the double-peak structure disappears when the size is reduced to 2 fm. This is similar to our results, that is, as q increases, the two peaks begin to merge into a broad one. And this phenomenon is very similar to c v /T 3 because the transport coefficient ζ is related to c 2 s and therefore also related to c v /T 3 . Besides, its response to q has a maximum near the pseudo-critical temperature T c .
Finally, the value of T 0 cannot be completely determined, so we consider T 0 as a free parameter to test the T 0 -dependence of our results. We find that T 0 has only a quantitative effect on the results, mainly manifested in that the pseudo-critical temperature T χ and T d decreases as T 0 decreases, as pointed out in Ref. [46]. However, our results are qualitatively independent of T 0 . Taking ǫ/T 4 as an example, as shown in Figs. 19, 20, we find that its change with q is not affected by T 0 .

IV. SUMMARY AND CONCLUSION
In this paper, combined with the Tsallis statistics and the PNJL model, we investigated the sensitivity of phase transitions, thermodynamic quantities, and transport coefficients to deviations from usual BG statistics. It was found that the chiral and deconfinement transition are still a crossover at finite temperature and zero quark chemical potential, independent of the non-extensivity parameter q. However, their corresponding pseudocritical temperatures T χ and T d decrease as q increases.
Regarding the influence of the parameter q on the thermodynamic quantities. On the one hand, we found that for ǫ T 4 , p T 4 and s T 3 , their high temperature limit is no longer the SB limit but the q-related Tsallis limit. But for c 2 s , due to a surprising cancellation, its high temperature limit is not affected by the parameter q. It should be noted that the finite-size effect does not change the SB limit. On the other hand, we found that as the q increases, the criticality of c v /T 3 and c 2 s will gradually or even completely disappear. This is consistent with the finite-size effect, where the same phenomenon occurs as the size decreases.
Under a constant value of relaxation time τ = 1/Γ = 1 fm, we calculated the transport coefficient as a function of temperature and found that η and σ rise monotonically with temperature increase, while ζ shows a maximum near the pseudo-critical temperature T c . About the influence of the parameter q on them, we found that η and σ increase as q increases at any fixed temperature, while ζ changes with q, similar to c v /T 3 . We also introduced a new susceptibility in order to study the response of thermodynamic quantities and transport coefficients to q in more detail. We found that their response patterns to q are different. For example, ǫ/T 4 , s/T 3 , ∆/T 4 and c 2 s have the highest response to q near the pseudo-critical temperature T c , while σ shows a maximum and a minimum values near the T c and then rises monotonically with temperature increase. For p/T 4 , its response increases with temperature increase and gradually stabilizes. Besides, it should be noted that the double-peak structure appearing in c v /T 3 and ζ is unique to Polyakov-loop potential U L . And as q increases, it gradually disappears. Interestingly, in the study of finite-size effects, the double-peak structure also disappears as the size decreases.
As a first step, we are only concerned with the zero quark chemical potential and finite temperature region and q > 1. Next, we will study the whole qdependence of the equation of state at finite temperature and finite quark chemical potential to further study its influence on the properties of protoneutron stars [65]. Furthermore, determining the existence and location of the critical end point (CEP) in QCD phase transition has been one of the main goals of relativistic heavy-ion collision experiments. For this purpose, the second phase of the beam energy scan at RHIC will be performed between 2019 and 2021 [66]. Therefore, the impact of non-extensive effect on the location of CEP is a question worthy of further study [4]. Moreover, Ref. [67] proposes an improved Polyakov-loop potential due to the backreaction of the quarks. Thus, it is also an interesting question to use the improved Polyakov-loop potential to study the non-extensive effect to compare with the usual Polyakov-loop potential. Finally, in order to qualitatively understand the influence of the parameter q on the transport coefficients, we have taken constant value of relaxation time in this present work. However, involved calculations of relaxation time at finite temperature and finite quark chemical potential in Tsallis statistics, incorporating different interaction channels might lead us to more realistic scenario. These issues are our future research directions.