Nonextensive Nambu Jona-Lasinio model of QCD matter revisited

We present a revisited version of the nonextensive QCD-based Nambu - Jona-Lasinio (NJL) model describing the behavior of strongly interacting matter proposed by us some time ago. As before, it is based on the nonextensive generalization of the Boltzmann-Gibbs (BG) statistical mechanics used in the NJL model to its nonextensive version based on Tsallis statistics, but this time it fulfils the basic requirements of thermodynamical consistency. Different ways in which this can be done, connected with different possible choices of the form of the corresponding nonextensive entropies, are presented and discussed in detail. The corresponding results are compared, discussed and confronted with previous findings.


Introduction
Some time ago we presented a nonextensive version of the QCD-based Nambu -Jona-Lasinio (NJL) model of manybody mean field theory describing the behavior of strongly interacting matter [1,2,3,4,5,6,7,8,9], the q-NJL model [10,11,12]. It was based on the nonextensive generalization of Boltzmann-Gibbs (BG) statistical mechanics used in the NJL model to its nonextensive version and was characterized by a dimensionless nonextensivity parameter q. At the same time the other many-body mean field theory of nuclear matter, the Walecka model [13,14,15] formulated on the hadronic level, has also been generalized to its nonextensive version [16,17,18], similarly as a number of other nonextensive approaches to different aspects of dense hadronic and quark matter, cf., for example, [19,20,21,22,23,24,25].
In short, introduction of nonextensivity to a given model consists in the replacement of the usual extensive Boltzmann-Gibbs (BG) distribution by a nonextensive Tsallis distribution [26,27,28,29] (preserving the original dynamical structure of the model 1 ): where exp q − X T a e-mail: jacek.rozynek@ncbj.gov.pl b e-mail: grzegorz.wilk@ncbj.gov.pl 1 As was done, for example, in the nonextensive hydrodynamical model discussed in [30,31].
X denotes the variable of interest, S is the BG-Shannon entropy and S q the Tsallis entropy, (for q → 1 both entropies coincide). However, one should be aware of the fact that the nonextensive distribution (1) can also be regarded as an example of a quasi-power law extrapolating between a pure, scale invariant, power-like form for large X and an exponential one for small X. In fact, it is observed in all branches of physics [26], including multiparticle production experiments [32,33,34,35,36,37,38,39,40,41,42,43,44,45], studies of complex systems [46,47,48,49,50] or systems on lattices [51,52,53,54]. As such, it can be derived in a variety of ways, not necessarily connected with statistical mechanics and not based on any entropy [55,56,57,58,59,60,46,47,48,49,50]. Because all our further considerations are based on nonextensive thermodynamics with some nonextensive form of entropy, it should be mentioned that such an approach is fully compatible with the usual traditional extensive thermodynamics [61,62,63,64,65,66].
Our motivation behind investigating the possible effects of introducing a nonextensive environment (with q not equal to unity) into the original extensive mean field approach of the NJL model remains the same as before [10,11,12]. Referring there for details, basically it is a desire to account, in a phenomenological way, for the numerous factors which make the assumptions of the mean field approach questionable. For example, it does not account for possible intrinsic (nonstatistical) fluctuations and correlations (especially short range ones), which are common when quark matter is produced in high energy heavy ion collisions. In such collisions it emerges in small and rapidly evolving samples, the spatial configuration of which is far from being uniform, and no global equilibrium is established [67,68,69,70]. The NJL model is therefore usually supplemented by some additional dynamical ingredients (like, for example, the Polyakov loops (cf., [71,72,73]) or by some retardation effects [5,6,7]). Instead of these, we propose to change summarily the type of environment by allowing it to be nonextensive, with the degree of nonextensivity provided by q − 1 and without specifying its dynamical origins. On the other hand, it is expected that addition to such an approach any dynamical ingredient of the type mentioned above should result in a diminishing of the value of |q − 1| used. In fact, a fitted value of q = 1 would signal that there is no need for any further nonextensivity based on nonextensive statistical mechanics 2 . Its description in terms of nonextensive statistical mechanics was replaced by some direct, dynamical one.
Unlike the simple situations leading to Eq. (1), in the q-NJL model one deals with fermions and has to account for both particles and antiparticles with nonzero chemical potentials. In such cases the formulation of the nonextensive formalism is not unique, especially when one wishes to fulfil the basic requirements of thermodynamical consistency [74]. First of all, the proper number of constituents in this case is not given by the usual nonextensive particle occupation numbers, n q , but rather by their q powers, n q q [75,76,77,78,79,80]. Further, the nonextensivity parameter q is not necessarily constant (actually this was the point overlooked in our previous work [10,11,12]). As shown in [74] the parameter q depends on whether one describes particles or antiparticles, (changing, for example, from q for particles to 2 − q for antiparticles), or, in other cases, it can depend on the value of the density n q (by changing, again, from q to 2 − q, but this time at the Fermi level at which n q = 1/2).
The effective number of constituents n q q was already used in [19,20,21,22,23,81,82] (with different choices of the form of nonextensive densities) and also in a recent new version of the q-Walecka model proposed in [83]. However, the necessity of using in some cases for the exponent both q and 2−q was fully explored in the analysis of dense hadronic matter only recently [84,85,86,87,88,89,90,91]. In this respect our previous calculations [10,11,12], performed assuming that the exponent q is always constant, were not correct. Therefore, we shall repeat them here, this time for all plausible choices of densities n q and for the correct choices of the corresponding n q q as proposed recently in [74]. At the same time the physical meaning of using the q > 1 or q < 1 type of nonextensive environment (usually connected with, respectively, intrinsic dynamical fluctuations [92,93,94] or correlations [95,96,97]) will also be discussed.
2 Such a situation was encountered some time ago in multiparticle production processes [34]. The gradual accounting for the intrinsic dynamical fluctuations in the hadronizing system by switching from a simple nonextensive distribution of the type of Eq. (1) to some dynamical formula with temperature fluctuations, substantially lowered the value of |q − 1|.
The outline of our work is as follows. After providing in Section 2 a short reminder of the NJL model used, we present in Section 3 the revisited version of our q-NJL model. Different ways of implementing it will be presented here and discussed in detail. Section 4 contains our results with special emphasis on the presentation and detailed discussion of the different forms of nonextensive entropy used and on their interrelations. In particular, we show that the usual ordering of entropies with respect to the value of the nonextensivity parameter q, S(q < 1) > S(q = 1) > S(q > 1) [26,27,28,29], in the case of q-NJL, with qdependent occupation numbers, can be reversed. We shall also present here results on the susceptibilities and the heat capacity of nonextensive nuclear matter. Section 5 summarizes our presentation.

Reminder of the NJL model used
With regard to the NJL model used by us, we follow the version presented in [8], with an effective lagrangian suitable for the bosonization procedure and with a four quark interaction only given by (for more details see [10]): wherem = diag (m u , m d , m s ) and It is invariant under the chiral SU L (3)⊗SU R (3) transformations described by coupling constant g S (except for the current quarks mass term) and contains a term breaking the U A (1) symmetry (described by thecoupling constant g D ), which reflects the axial anomaly in QCD. D abc are the SU(3) structure constants d abc for a, b, c = (1, 2, . . . , 8) whereas D 0ab = −δ ab / √ 6 and D 000 = 2/3. We work with q = (u, d, s) quark fields with three flavors, N f = 3, and three colors, N c = 3, λ a are the Gell-Mann matrices, a = 0, 1, . . . , 8 and λ 0 = 2/3 I.
Integrating over the momenta of quark fields in the functional integral with L ef f one obtains an effective action expressed in terms of σ and ϕ, the natural degrees of freedom of low energy QCD in the mesonic sector (Tr stands for taking the trace over indices N f and N c ): The first variation of W ef f results in the gap equations for the constituent quark masses M i : with cyclic permutation of i, j, k = u, d, s and with the quark condensates (S i (p) is the quark Green function, m i is the current mass of quark of flavor i; notice that nonzero g D introduces mixing between different flavors). We consider a system of volume V , temperature T and i th quark chemical potential µ i characterized by the baryonic thermodynamic potential of the grand canonical ensemble, with quark density equal to ρ i = N i /V , baryonic chemical potential µ B = 1 3 (µ u + µ d + µ s ) and baryonic matter density ρ B = 1 3 (ρ u + ρ d + ρ s ). The internal energy, E, the entropy, S, and the particle number, N i , are given by [8,9] (here E i = M 2 i + p 2 ): All these quantities depend on the quark and antiquark occupation numbers, n i andn i . They can be obtained in different ways; we shall use the standard Jayne's information-theoretic approach and derive them from a given measure of information represented by the entropy functional (13). Extremalizing it under constraints imposed by the total number of particles,N , and total energy of the system,Ê (E i is the energy of the i-th energy level) [98,99,75], i (n i −n i ) =N and i (n i +n i ) E i =Ê, (15) one finds that (e(x) = exp(x)): The values of the quark condensates in Eq. (8) are now given by Eqs. (8) and (18) form a self consistent set of equations to calculate, for a given T and µ, the effective quark masses M i and values of the corresponding quark condensates, q i q i , and, further, all quantities of interest. With these effective masses M i the occupation numbers (16), the form of which was obtained assuming a noninteracting Fermi gas, will get their dynamical input defined by the form of lagrangian used, Eqs. (4) and (7). The model is fixed by the coupling constants g S and g D and by the cutoff in three-momentum space Λ, which is used to regularize the momentum space integrals and the current quark masses m i . The values of all parameters are the same as in [10].
3 The q-NJL model revisited

The problem of thermodynamic consistency
Thermodynamical consistency means that one has to preserve the standard thermodynamical relationships among thermodynamical variables such as, for example, entropy, energy and temperature, As shown in [100,101], any thermostatistical formalism constructed by following Jayne's maximum entropy prescription complies with the thermodynamical relationships. A thermodynamically consistent formulation is then obtained by the appropriate identification of the relevant constraints with the extensive thermodynamical quantities (like number of particles N or energy E) and by identification of the corresponding Lagrange multipliers with the appropriate intensive thermodynamical quantities (like temperature T and chemical potential µ). In our case we shall extremalize some appropriately chosen entropic measures using some specifically chosen constraints. Because this procedure is not unique, we shall discuss three different versions of the q-NJL model, based on different choices of entropic measureS q proposed in [74].
3.2 Approach with Tsallis' cut-off prescription for q > 1 For q > 1 the straightforward approach uses a q-generalization of the entropic measureS in Eq. (13), which reduces to it for q → 1, where ln q x = To fulfil the basic requirements of thermodynamical consistency with such a form of entropy functional, the constraints must have the following form [74]: Such a form reflects the fact that in the nonextensive environment parameterized by q the relevant number of constituents is given not by n q , but by n q q . This is because, as shown in [46,47,48,49,50], the occurrence of the nonextensivity can be formally understood as a result of the action of some effective interactions between constituents of the system, which therefore lose their independence. The strength of these interactions is proportional to ζ = q − 1: they are attractive for ζ < 0 and repulsive for ζ > 0, and they are not accounted for by the usual description of the system considered; for example, in our case of the q-NJL model they act on top of the interaction described by the lagrangian defined by Eq. (4).

Extremalization ofS
(a) q under constraints (22) results in the following quark and antiquark occupation numbers 3 : and As in the extensive case, Eq. (16), the above effective occupation numbers are derived for a Fermi gas, i.e., for a situation with no interactions. However, this time it is immersed in a nonextensive environment, which, as mentioned above, introduces some effective interaction [46,47,48,49,50]. Therefore, it is in reality an assumption that we shall use them in further calculations where E qi = M qi + p 2 and M qi denotes the quark mass obtained from the q-version of the gap equation (8), cf. Eq. (52) below. In addition to the dynamical input mentioned before when discussing Eq. (16), this is additional source of nonextensivity in a q-NJL model. For q → 1 we recover Eqs. (16), (8) and (18) of the usual NJL model. The same remark applies to the other two approaches presented below 4 .
Definition (25) must be supplemented by a condition ensuring that thefunction e q (x) is always non-negative real valued (known as Tsallis' cut-off prescription), There are no limitations on the distribution of antiparticles,n qi . We observe always that n q>1 (x) > n q=1 (x) (cf. also [76,77]). 3 In what follows we shall use the notation nq or nq (and eq or eq) in situations where q in the above definitions will be replaced by a suitably definedq orq. In other cases the q in Xq are just subscripts indicating the nonextensive origin of variable X. 4 Note that for T → 0 one always gets nq(µ, T ) → n(µ, T ), irrespective of the value of q. This means that we can expect a nonextensive signature only for high enough temperatures. A similar situation is encountered in models using a Polyakovloop potential (for example, in a Polyakov-loop extended quark meson model [102]).

Approach with Tsallis' cut-off prescription for q < 1
The case where q > 1 is the usual range explored in all previous investigations of this type, cf. [16,17,18,19,20,21,22,23,81,82]). However, because in the q-NJL model we must consider on the same footing both particles and antiparticles, we have to address how to formulate our approach for the q < 1 case. It is easy to check that using the previous description we would have to replace in this case the condition (26) by . (30) However, this means that for reasonable temperatures T the phase space for antiparticles in such an approach would be practically completely cut-off. Therefore, if we intend to extend this approach to the q-NJL model and consider the case q < 1, we have to modify the form of the entropy functional and replaceS (a) by [74] with accordingly modified constraints, The resulting occupation numbers are still given by Eq. (23) with nonextensive parameter q for particles and with q replaced by its dual,q = 2 − q, for antiparticles. This means that innq i one now has In this case we always observe that n q<1 (x) < n q=1 (x). With such a choice of nonextensive parameter there are no longer any problems for antiparticles with regard to phase space limitations and their physical meaning. However, the price one pays is that now particles and antiparticles are described by different, albeit connected, nonextensive parameters: q for particles and its dual,q = 2 − q, for antiparticles. Therefore the choice q < 1 for particles meansq > 1 for antiparticles (i.e., it remains the same as in Section 3.2 above). The possible justification of such a choice could be that antiparticles, in contrast to particles, do not have a Fermi level (in the sense that their states are not fully occupied). Therefore, they are not so strongly correlated as particles, but rather experience dynamical fluctuations (believed to be described by q > 1 [92,93,94]).

Approach without Tsallis' cut-off prescription
The Tsallis cut-off prescriptions needed in the above two choices of the entropy functional,S (a) q in Eqs. (20) and S (b) q in Eq. (31), limit considerably the allowed phase space and, concerning the q-NJL model, they apparently do not have any physical justification. Therefore, in [76,77,78] an alternative approach was proposed eliminating the cut-off prescription by a suitable choice of the entropy functional. Thus far it has been used in [79,80,83,84] for particles only. In the case of the q-NJL model where both particles and antiparticles have to be considered together this approach has to be modified. The proposed choice can be summarized as follows [74]: -for fermions above the Fermi sea level, i.e. for x qi > 0 (or n qi ∈ [0, 1/2]), we always use q > 1; -for fermions below the Fermi sea level, i.e. for x qi < 0 (or n qi ∈ 1 2 , 1 ), we always useq = 2 − q < 1; -for antifermions we assume q > 1 over the whole region of x (or for n qi ∈ [0, 1]). These requirements correspond to the following, third, choice of the corresponding q-entropy functional: and constraints in the following form: (35) whereq now changes at the Fermi surface in the following way:q With such choices one has the following occupation numbers: for particles : for antiparticles :n qi = 1 e q (x qi ) + 1 (38) where In this case, instead of being either always greater (for q > 1, as in case (a)) or smaller (for q < 1, as in case (b)) than the respective values for the extensive case, we observe that nq < n q=1 below the Fermi surface (i.e., for x < 0) and nq > n q=1 above the Fermi surface (cf. also [76,77]. The most important observation which must be kept in mind is that whereas in this prescription nq i has smooth behavior for all values of x qi , theq (36) and therefore also nq qi have a jump at x qi = 0, i.e., for n qi = 1/2 (at Fermi surface), the physical meaning of which is so far unclear.

Problem of particle-hole symmetry in a nonextensive environment
In an extensive situation (at least in the mean field approximation), with occupation numbers given by Eq. (16), the particle-hole symmetry is preserved, and one can interpret holes among the negative energy states as anti-particles with the corresponding positive energy. However, any short range interaction violates relation (39) by enhancing the tail of the distribution above the Fermi level [103]. In a nonextensive environment this symmetry (39) is violated and holds only in a dual fashion [104]. For example, when implementing nonextensivity as in Sections 3.2 and 3.3, we observe that On the other hand, the prescription described in Section 3.4 apparently restores the particle-hole symmetry, because now However, one should keep in mind that, as was stressed before, the relevant momentum distributions of particles are now given by n q q , not by n q , because one deals not with free particles but with particles embedded in nonextensive environment [46,47,48,49,50]. It means therefore that in our case this symmetry does not work.

Basic formulas of the q-NJL model
The basic quantity is now the q-version of the grand canonical potential (10), The pressure and the energy density are defined as, respectively, where Ω q (0, 0) = E q (0, 0) denotes the vacuum energy. The first derivatives give, respectively, q-entropy and q-density, The second derivatives result in nonextensive versions of the heat capacity, C µ , and the barionic susceptibility χ B , Because one can define temperature the T and the heat capacity C µ as, respectively, The corresponding nonextensive energy, E q , entropy, S q , and number density, N q , are now given by N qi = N c π 2 V p 2 dp n q qi −n q qi .
The q-version of the gap equation (8) for the constituent quark mass in the nonextensive environment is now with a q-dependent quark condensate, which, following [19,20,21,22,23,81,82], is adopted in the form: Depending now on the form of the entropic functional S (R) q used, there are the following three possibilities for the nonextensivity parameters.
The first choice is that for q > 1 as given by Eq. (20) supplemented by constraints (22) in which the parameter q is constant and remains the same for particles and antiparticles. As a result, the n qi andn qi are given by Eq. (23), and the effective distributions, n q qi andn q qi , have the same value of the power index q for all values of n qi andn qi .
The second choice, to be used if one intends to investigate the case of q < 1, is as given by Eq.(31) with constrains (32). In this case the parameter q for particles is replaced for antiparticles by its dual,q = 2 − q. As a result, whereas formally n qi andn qi are again given by Eq. (23), the parameter q for particles is replaced byq for antiparticles. The same situation is encountered for the effective distributions which are now equal to, respectively, n q qi andnq qi . The third choice, for q changing at the Fermi surface, is thatS as given by Eqs. (34) with constraints (35) and with n qi andn qi given by Eqs. (37) and (38). The effective occupation numbers are now nq qi (withq defined by Eq. (36)) andn q qi . Note that this time one encounters for particles a discontinuity in this variable on the Fermi surface, i.e. for n qi = 1/2. However, there is no such discontinuity for antiparticles.

Introductory remarks
We first check numerically the thermodynamical consistency of q-NJL. We have the following set of differential equations: where the energy density ε q = E q /V , the entropy density s q = S q /V and the density ρ q = N q /V . The relations to be checked are We have checked these relations taking as input µ = 322 MeV and T = 60 MeV, for case (a) with q = 1.1, case (b) with q = 0.9 and case (c) with q = 1.1, and also for the extensive case q = 1. We have also checked that these relations are satisfied even for a broader range of values of the parameter q, ranging from q = 0.8 to q = 1.2. They are all satisfied with an accuracy of better than 1%. Our previous formulation of the q-NJL model [10] followed essentially choice (c) with regard to calculations of the nonextensive occupation numbers n q andn q , but their effective values, n q q andn q q , were wrongly used with a constant value of the exponent q. The three formulations of the q-NJL model presented in this work, withS In what follows we provide such a comparison, calculating the dependencies on chemical potential µ and temperature T of the entropy S q , pressure P q , normalized density ρ/ρ 0 , heat capacity C µ and barionic susceptibility χ B .

Results for entropiesS (R=a,b,c) q
We start with a presentation of the entropiesS (R) q corresponding to different ways (R = a, b, c) of introducing nonextensivity into the q-NJL model. They are shown as functions of the scaled variable xT /µ in Figs. 1 and 2. The results shown will allow for a better understanding of our further results in Section 4.3. In Fig. 1 we present a schematic view ofS multiplied by square momentum, p 2 , and with running mass M qi , will be further integrated and result in entropy S q in Eq. (50), which will be shown below in Fig. 3.
Some comments are in order here. It is known that in the case when occupation numbers do not depend on q one observes the ordering of entropies from subextensive for q > 1 to superextensive for q < 1 withS q>1 <S q=1 < S q<1 [26,27,28,29]. However, it turns out that in the case where occupation numbers depend on q, like in our q-NJL model, the situation is more complicated. Two features can be observed. First, this ordering changes with x. It starts withS q , but somewhere above the Fermi surface it reaches the point (marked by the vertical line) where all these entropies coincide and, from there, their order is reversed. For M u = 0 this point is located at x = 0.25, but as seen in Fig. 2 it shifts slowly to higher values of x with increasing mass M qi . At the same time the minimum value of x (corresponding to p = 0 and equal to M 2 u + p 2 1/2 − µ /T , cf. Eq. (24)) shifts with increasing M u towards higher values, nearer the point where all entropies coincide and their order is reversed. This effect is strongly enhanced for higher momenta and, eventually, almost all contributions to the final entropy will come from this region. This will result in the ordering of entropies seen below in Fig. 3. The same shift would be observed for case (c) (not shown here since, as was stated previously, the results for this case can be obtained by a suitable composition of the results for cases (a) and (b)). This means that for higher masses and momenta the dominant component will therefore be given by q (a) = 1.1 (which will show up in Fig. 3 below).
We close this section by noting that it can be checked that such behavior of the entropiesS (R) q is caused by differences in the tails of the occupation number distributions n q (x) as functions of x, additionally enhanced by the fact that in the definition of differentS (R) q they enter as powers of q, n q q (x). For x = 0 all distributions coincide to n q (0) = 1/2. However, in case (c), where at x = 0 the power index q changes to 2 − q (for example from q = 1.1 to 2 − q = 0.9) n q q (0) = n (2−q) q (0). All these results were calculated neglecting the contribution of the sea (i.e. neglecting antiparticles). We have checked that their contribution does not change the conclusions reached here, their effect being of the order of 1% only.

Results for entropy, pressure and compression
After the preparatory explanations in the previous section we present now results for entropy S (Fig. 3), pressure P (Fig. 5) and density ρ (Fig. 4). Calculations were performed for T = 60 MeV and µ = 322 MeV, for different values of q corresponding to different realizations of the q-NJL model and compared with the BG case of q = 1. This choice of values of T and µ is dictated by the fact that, as can be seen in Table 1, they are near to the critical values corresponding to the values of nonextensive parameter q considered here. For this choice the critical effects displayed in the presented figures are most visible. Chemical potentials and temperatures at which one observes rapid changes of entropy (Fig. 3) and small jumps in the pressure (Fig. 5), correspond to regions of phase transition in which the density of fermions changes rapidly, approximately from ρ 0 to 2ρ 0 , as seen in Fig. 4.
The general characteristic feature of these results is the horizontal shift of the phase transition point and of all curves for entropy S and density ρ/ρ 0 , when calculated for a fixed value of temperature T = 60 MeV, towards smaller values of µ for q < 1 (case (b)) and towards greater values of µ for q (a) = 1.1 (case(a)). For densities below and above well separated critical points all curves converge. At the same time S Behavior of pressure P is shown in Fig. 5. Its increase as a function of µ is dictated by the derivative ρ = ∂P/∂µ| T  (and depends on the behavior of ρ presented in Figure 4), whereas its increase as a function of T is given by the derivative s = ∂P/∂T | µ (and is given by the behavior of S in Fig. 3). Looking at Figs. 3-4 from this perspective, we observe generally that both the entropy and density are lower below the phase transition than after it. The corresponding increase of pressure before the phase transition is therefore weaker than in unconfined phase. This behavior is more pronounced when considered as a function of temperature T rather than as a function of chemical potential µ. Note that because, as seen in Fig. 4, densities are very similar for different realizations of the q-NJL model, the resulting pressure curves as a function of chemical potential are almost parallel. Because entropies differ much more strongly the dependence on temperature is more divergent for different realizations of the q-NJL model.

In the vicinity of the phase transition point
We shall present now the behavior in the vicinity of critical values of T and µ and present results for pressure, compression, heat capacity C µ , baryon number susceptibility χ and the phase diagram. Fig. 6 presents the pressure at the critical temperature T cr (critical isotherms) for the corresponding critical values of the chemical potential, µ cr , as a function of compression ρ/ρ 0 . The values of T cr and µ cr determine the critical values of pressure and density corresponding to the inflection points for which the first derivative of pressure by compression vanishes. They are listed in Table 1. Remarkably, in all these cases the critical density remains practically the same. The other interesting observation is that the results for choice (a) essentially coincide with those for choice (c). For q > 1, we observe an increase of the pressure P (corresponding to a similar increase of the entropy in this case observed in Fig. 3). This increase may simulate some additional repulsion between massive quarks constituting nucleons at        short distances (which is present in quark-meson coupling models based on the mean field approach [106,107]). On the contrary, for q < 1 we observe a decrease of the critical pressure (corresponding to a decrease of entropy, seen in Fig. 3). This could simulate confinement of quarks in nucleons which introduces restrictions in phase space by changing their Fermi motion.
In Fig. 7 we present compression as a function of chemical potential µ and temperature T but now, unlike in Fig.  4, calculated in the vicinity of the phase transition point (the corresponding values of T cr and µ cr are given in Table 1). Because both T and µ now change their values, the observed behavior of ρ/ρ 0 is very different from the previously one. It is dictated by the fact that, as seen in Table 1, in our q-NJL model we have T  In Fig. 8 we present our results for, respectively, heat capacity (not calculated before in [10,11,12]) and baryon susceptibilities. Both were calculated for the critical values of the chemical potential and temperature listed in Table  66 67  In this region we also observe that the chemical potential at the critical point is lower, µ crit = 312 MeV, in comparison to its extensive value for q = 1, µ crit = 318.4 MeV. As a result, the critical points for q (b) < 1 lie very close the extensive phase border, as seen in Fig. 9. For q (a) = 1.1 (and partially also for q (c) = 1.1) the critical temperature T crit remains very near its extensive value for q = 1, T crit = 70.5 MeV, but the critical chemical potential is shifted towards larger values, see Fig. 9 and Fig. 8 (right panel). As seen in Fig. 7 and Fig. 9 the results for T crit are for choice (c) are very near to those for choice (a) whereas the results for µ crit are for choice (c) very near to the results for q = 1. Consequently, whereas the critical coordinates in cases (a) and (b) depart from the extensive phase border line in Fig. 9, the critical coordinates for choice (c) are very close to the extensive limit of q = 1. This reflects the fact that in choice (c) we use at the same time both q > 1 and q < 1 (depending on the region of phase space actually considered) minimizing therefore the departure from the extensive dynamics.
The results presented in Figs. 9, 6 and in the right panel of 8 can be compared with the corresponding results obtained in our previous version of the q-NJL model (with Figs. 4 and 8 in [10], with Fig. 1 in our first paper in [11,12] and with Figs. 1-3 and 5 in our second paper there). This allows us to estimate what kind and how big are the changes introduced by proper treatment of thermo-dynamical consistency. Quantitatively the corresponding results look similar, there are differences in what concerns sensitivity to deviation from extensivity leading to quantitatively similar effects: the previous |q − 1| ∼ 0.02 should be confronted with the present ∼ 0.1. The detailed shape of the phase diagram or detailed shapes of some distributions are also different, although, not dramatically so. The most visible, but still rather small, difference concerns the case with q < 1. For example, the previously presented χ was much smaller and broader than the corresponding result for q > 1, at present they are comparable in shape and height (compare right panel of our Fig.8 and Fig.5 in the second paper of [11,12]). Also, the sensitivity to changes in q in the vicinity of the phase transition point seems to be reversed. Comparing our Fig. 6 and Fig. 1 in [11,12] one notices greater changes for q < 1 than for q > 1 observed before, and the reverse behavior is seen in the present results. This is so far the most visible difference resulting from the fact that now q < 1 is described entirely by case (b) only. The most similar to what was used before in [10,11,12] is case (c), already explained in detail.

Summary and conclusions
The results presented above confirm the previously reached conclusions concerning the applicability and usefulness of the nonextensive approach as one of the possible extensions of the original mean field theory, extending considerably the range of its applicability. However, this time they were obtained in compliance with all the conditions that must be satisfied when implementing a nonextensive approach to systems with particles and antiparticles and a nonzero chemical potential (which were so far overlooked) [74]. Out of three possible formulations investigated in detail in sections (3.2), (3.3) and (3.4) (and denoted by (a), (b) and (c), respectively), the most consistent, in our opinion, is approach (a) for q > 1 or (b) for q < 1.
Approach (c) requires more attention. There are two features which should be commented on. The first is that our results in Fig. 3 show that the jump in entropy S (c) (for q (c) = 1.1) occurs almost at the same point as that for the extensive case of q = 1. This is because in case (c) we have, in fact, two values of the nonextensive parameter: q above the Fermi surface and 2 − q below it (equal, correspondingly, to 1.1 and 0.9 in our case). As can be seen in Fig. 3, the corresponding shifts for q > 1 and q < 1 are in opposite directions. Therefore, we see here a kind of cancellation of two opposite effects. The second feature, visible in Fig. 1 (lower panel) is the jump in entropy on the Fermi surface (which occurs because of the above mentioned change in nonextensivity parameter there). However, in nuclear matter we do not observe a sharp change of entropy on the Fermi surface (i.e. for x = 0) in the extensive situation of q = 1 which could possibly be modeled by introducing a nonextensive envi-ronment with q = 1. The transition from the bound state to the continuum is always gradual 6 [110,111].
As discussed before in Section 3.4, case (c) was proposed in order to avoid Tsallis' cut-offs introduced in cases (a) and (b) (cf. Eq. (27)). These cut-offs result in the artificial fixing of the respective occupation numbers in some regions of phase space to be equal either to unity (particles) or zero (antiparticles) (cf. Eqs. (26) and (27) or (28)-(30)). They are therefore regarded as unjustified additional ingredients of the q-NJL model. However, it seems that the price of introducing case (c), which eliminates the need for such restrictions, is too high to be acceptable for the reasons mentioned above. The cut-off of part of the phase space seems to be more reasonable in this respect. To add to these arguments, note that when crossing the Fermi level the quark mass does not change, whereas it changes in the phase transition.
In our revisited non-extensive, but thermodynamically consistent, q-NJL model, we can perform calculations for both q < 1 and q > 1 statistical effects. Some comments regarding their possible roles are therefore necessary. Note first that the main approximation in the NJL mean field theory description concerns the phase with condensates, where the quarks are not confined but rather become less massive approaching the phase transition region and practically massless when crossing it (in the same way as in our previous version of the q-NJL model, cf. Fig. 1 in [10]). In the usual NJL mean field model confinement is not present; it can be introduced by adding a dynamical gluon field, for example, in the form of a Polyakov loop [71,72,73]. On the other hand, confinement of quarks in nucleons introduces restrictions in the allowed phase space by changing their Fermi motion and decreasing the chemical potential. It increases quark arrangement in nuclear matter and, consequently, decreases the entropy and effectively also the pressure between such composite objects as nucleons. Such effects are visible in our q-NJL model with q < 1. This means that for the phase with vacuum condensates the q < 1 description could adequately describe the correlations which are responsible for the changes in the arrangement of quarks and therefore model, to some extent, the effect of confinement. These correlations will survive even in the unconfined phase (i.e. entropy there remains lower than in other cases). Also, as seen in Table 1, the critical temperature for q < 1 is higher (T (b) cr = 70.5 MeV), whereas the corresponding entropy becomes smaller. This is understandable because at 6 A singularity at the Fermi surface reported recently in [108] occurs only as a quantum effect at the vanishing temperature. On the other hand, such behavior is observed in solid state physics, cf. for example [109]. In nuclear matter one could think of replacing the jump in q in approach (c) from q to 2 − q by some smooth transition taking place between, say, x − ǫ and x + ǫ. In fact, similar propositions in this direction, using a suitable modifications of the occupation numbers for x < 0, were already presented and discussed in [84,104]. This, however, requires a completely new and demanding analysis of the nonextensive approach, which is beyond the scope of the present work. lower entropy we need larger temperature to convey the same amount of energy.
For the case of q > 1 (case (a)) the average single particle energy (chemical potential) is shifted towards higher values and the entropy and pressure are also higher, increasing the Fermi energy. The critical temperature (both for case (a) and case (c)) remains practically the same as for the extensive case with q = 1. The increase of pressure may simulate some additional repulsion (at short distances) between massive quarks constituting nucleons (see, for example, the quark-meson coupling models [106,107]). This means that, effectively, q > 1 could be regarded as emerging from the increasing nonstatistical fluctuations found in the decay of dense fireballs produced in heavy ion collisions and in other multiparticle production experiments [32,33,34,35,36,37,38,39,40,41,42,43,44,45].
Summarizing, the nonextensive description accounts (in a phenomenological way) for all situations in which one expects some dynamical correlations in the quarkantiquark system considered in our q-NJL model. It must be stressed that the nonextensive approach is not a substitute for any part of the interaction described already by the lagrangian of the NJL model, cf. Eq. (4). It rather provides a different environment which can have some dynamical effects, so far undisclosed but simply parameterized by nonextensivity q. From this perspective case (b) with q < 1, which has lower entropy, seems to be more suitable to describe the nonextensive mechanism in the Equation of State (EoS) of dense hadronic matter in a restricted phase space, like, for example, in (proto)neutron stars. Respectively, the case (a), with q > 1 and with higher entropy, is more suitable for situations in which one expects dynamical fluctuations, for example to describe heavy ion collision where dense nuclear matter is created out of equilibrium and quickly decays into hadrons 7 . Our work therefore presents arguments that the critical properties of nuclear matter in two different environments can be different, although the phase transition occurs at the same density or compression. The critical temperature is higher for nuclear matter created in (proto)neutron stars but the critical value of the chemical potential will be bigger for nuclear matter created in heavy ion collisions 8 We conclude with two remarks. First: our approach should not be confounded with the similar in spirit approach based on quantum algebras (or on the so called q-deformed algebras) which was used to formulate a qdeformed NJL model [114]. Their common feature is the use of some suitable deformation of the mean field NJL model (based on nonextensive statistical mechanics in our case and on quantum algebras in [114]), which may ac-7 This fully agrees with the expected meaning of the parameter q discussed previously (in [95,96,97] for q < 1 and in [92,93,94] for q > 1). 8 The application of the formalism presented here to calculations of the EoS for neutron matter with some admixture of protons, including also hyperons for higher densities, and to estimate the upper limit for the mass of neutron stars [112,113] is underway. It is intended to continue recent investigations presented in [82,105]. count for intrinsic correlations and fluctuations that go beyond the mean field formulation and, in a certain limit, approach the more realistic lattice calculations. Second: there exists another, potentially very interesting, approach to nonextensivity, based on the so called Kaniadakis entropy [115,116,117]. In [118] it was used to study the formation of the quark-gluon plasma formation (and compared with nonextensive approach) whereas in [119] it was used to investigate the relativistic nuclear EoS in the context of the Walecka quantum hadrodynamics theory.