Nonextensive effects on QCD chiral phase diagram and baryon-number fluctuations within Polyakov-Nambu-Jona-Lasinio model

TL;DR

This paper introduces a nonextensive statistical mechanics version of the PNJL model to study QCD phase transitions and baryon fluctuations, revealing non-monotonic CEP behavior and aligning with experimental data.

Contribution

It develops a nonextensive PNJL model incorporating the parameter q, extending the traditional model to account for non-BG effects and analyzing their impact on the QCD phase diagram.

Findings

CEP moves non-monotonically with q

Results align with experimental data for q=1.07

Baryon fluctuation moments are consistent with observations

Abstract

In this paper, a version of the Polyakov-Nambu-Jona-Lasinio (PNJL) model based on nonextensive statistical mechanics is presented. This new statistics summarizes all possible factors that violate the assumptions of the Boltzmann-Gibbs (BG) statistics to a dimensionless nonextensivity parameter $q$, and when $q$ tends to 1, it returns to the BG case. Within the nonextensive PNJL model, we found that as $q$ increases, the location of the critical end point (CEP) exhibits non-monotonic behavior. That is, for $q<1.15$, CEP moves in the direction of lower temperature and larger quark chemical potential. But for $q>1.15$, CEP turns to move in the direction of lower temperature and lower quark chemical potential. In addition, we studied the moments of the net-baryon number distribution, that is, the variance ($\sigma^{2}$), skewness (S), and kurtosis ($\kappa$). Our results are generally consistent with the latest experimental data, especially for $\sqrt{S_{NN}}>19.6\ \mathrm{GeV}$, when $q$ is set to $1.07$.

Tables

Table 1. Table 1: Parameter set used in our work.

$a_0$	$a_1$	$a_2$	$b_3$	$T_0\left(\mathrm{MeV}\right)$
3.51	-2.47	15.2	-1.75	190

Table 2. Table 2: Parameter set used in our work.

$\Lambda \left(\mathrm{MeV}\right)$	$g\left({\mathrm{MeV}}^{-2}\right)$	$m\left(\mathrm{MeV}\right)$	${\alpha }_1$	${\alpha }_2$
631.5	$5.498\times {10}^{-6}$	5.5	0.2	0.2

Table 3. Table 3: Numerical results of CEP for different nonextensivity parameter q 𝑞 q . The unit of T c subscript 𝑇 𝑐 T_{c} and μ c subscript 𝜇 𝑐 \mu_{c} is MeV.

$q$	$1$	$1.05$	$1.08$	$1.1$	$1.12$	$1.15$	$1.17$	$1.2$
$T_c$	159.5	135	119.8	115.7	113	109.4	107.7	104.7
${\mu }_c$	171	251.4	284.9	288.4	290	290.2	287.4	281

Equations

$${\rho }_{BG}(x)=C\exp (x)⟶{\rho }_q(x)=C_q{\exp }_q(x),$$

$${\exp }_q(x)=[1+(1-q)x{]}^{\frac{1}{1-q}},$$

$${\ln }_q(x)=\frac{x^{1-q}-1}{1-q}.$$

$${\mathcal{L}}_{PNJL}$$

$$\Phi =\frac{⟨T{r}_{c}L⟩}{N_c},\bar{\Phi }=\frac{⟨T{r}_c{L}^{†}⟩}{N_c},$$

$$L(\vec{x})=\mathcal{P}\exp (i{\int }_0^{\beta }{A}_4(\vec{x},\tau )d\tau ).$$

$$\frac{{\mathcal{U}}_{\mathcal{L}}}{T^4}=-\frac{a(T)}{2}{\Phi }^2+b(T)\ln [1-6{\Phi }^2-3{\Phi }^4+8{\Phi }^3],$$

$$a(T)$$

$$b(T)$$

$$G=g[1-{\alpha }_1{\Phi }^2-2{\alpha }_2{\Phi }^3].$$

$$\Omega$$

$$\Omega (\mu ,T,M,\Phi )$$

$$M=m-2G(\Phi )\sigma ,$$

$$F^{\pm }$$

$$\frac{\partial \Omega }{\partial M}=\frac{\partial \Omega }{\partial \Phi }=0.$$

$${\chi }_T=\frac{\partial \sigma }{\partial T}.$$

$${\Omega }_q$$

$${\Omega }_q(\mu ,T,M,\Phi )$$

$$F_q^{\pm }$$

$$[1+(1-q)x]>0.$$

$$M$$

$$0$$

$$n_q(T,\mu )=\frac{e_q^q(-3(E_p-\mu )/T)+\Phi (e_q^q(-(E_p-\mu )/T)+2e_q^q(-2(E_p-\mu )/T))}{[1+3\Phi (e_q(-(E_p-\mu )/T)+e_q(-2(E_p-\mu )/T))+e_q(-3(E_p-\mu )/T){]}^q},$$

$${\bar{n}}_q(T,\mu )=\frac{e_q^q(-3(E_p+\mu )/T)+\Phi (e_q^q(-(E_p+\mu )/T)+2e_q^q(-2(E_p+\mu )/T))}{[1+3\Phi (e_q(-(E_p+\mu )/T)+e_q(-2(E_p+\mu )/T))+e_q(-3(E_p+\mu )/T){]}^q}.$$

$$S\sigma$$

$$\kappa {\sigma }^2$$

$${\chi }_B^{(n)}=\frac{{\partial }^{n-1}}{\partial {\mu }_B^{n-1}}{\rho }_B=\frac{{\partial }^{n-1}}{3^n\partial {\mu }^{n-1}}\rho (T,\mu ).$$

$$\rho (T,\mu )$$

$$T$$

$${\mu }_B$$

Full text

Key-words: nonextensive statistics, Polyakov-Nambu-Jona-Lasinio model, QCD phase diagram, baryon-number fluctuations.

PACS Number(s): 12.38.Mh, 12.39.-x, 25.75.Nq, 12.38.Aw

I INTRODUCTION

Since the discovery of quark-gluon plasma (QGP) at the Relativistic Heavy-Ion Collider (RHIC) and Large Hadron Collider (LHC), the determination of QCD phase diagram at high quark chemical potential, especially the search for the location of critical end point (CEP), has been the focus of experiments, such as the BES-II program at RHIC and the CBM experiment at FAIR. Theoretically, due to the “sign problem” encountered by lattice QCD, people use various effective models to study the QCD phase diagram, such as chiral perturbation theory Espriu et al. (2020), finite energy sum rules Ayala et al. (2011), Dyson-Schwinger Equations (DSEs) Fischer (2019); Xu et al. (2015); Zhao et al. (2019), Nambu-Jona-Lasinio (NJL) model and Polyakov-Nambu-Jona-Lasinio (PNJL) model Buballa (2005); Cui et al. (2018); Fuseau et al. (2020); Zhao (2020); Zhao et al. (2020).

It is worth noting that Boltzmann-Gibbs (BG) statistics are often used in the study of QCD phase diagram. However, in relativistic heavy ion collisions, this is not satisfactory. The collision system is small enough and evolving rapidly that global equilibrium is difficult to establish. As a result, the probability distribution is no longer the standard $e$-exponential distribution in BG statistics, but becomes the so-called power-law tailed distribution (i.e. the Tsallis distribution introduced next), which is supported by more and more relativistic heavy ion collision experiments Sharma and Kaur (2019); Bhattacharyya et al. (2018); Abelev and et al. (2007); Adare and et al. (2011); Aamodt and et al. (2011); Aad and et al. (2011); Khachatryan and et al. (2010).

The most notable difference between Tsallis statistics and BG statistics is that the former replaces the probability distribution as follows Osada and Wilk (2008); Rożynek and Wilk (2016); Ke-Ming Shen and Wang (2017):

[TABLE]

where

[TABLE]

and the inverse function is

[TABLE]

As $q\rightarrow 1$, BG statistics is recovered.

In fact, the Tsallis distribution appears not only in relativistic heavy ion collision systems, but also in complex systems with non-ergodicity, multifractality, long-range correlations, long-term memory, etc Tirnakli and Borges (2016); Cirto et al. (2018); Tsallis and Cirto (2013). Ref. Tirnakli and Borges (2016) shows in a very clear manner that due to a break in ergodicity, the system has crossed from BG statistics to Tsallis statistics. Furthermore, it should be noted that the research in Ref. Deppman et al. (2020) shows that the renormalizable field theories lead to fractal structures, which can be studied by using Tsallis statistics. Therefore, in this paper we will use Tsallis statistics to study the QCD phase diagram. And the paper is organized as follows: In Sec. II, we generalize the PNJL model to its nonextensive version. And in Sec. III, we focus on the impact of nonextensive effects on the QCD phase diagram, especially on the location of CEP. In addition, we study the moments of the net-baryon number distribution and compare them with latest experimental data. Finally, we give a brief summary of our work in Sec. IV.

II PNJL and nonextensive pnjl model

II.1 PNJL model

The Lagrangian of the two-flavor and three-color PNJL model for equal-mass quark is Ratti et al. (2006); Fukushima (2004):

[TABLE]

where $m$ is the current quark mass, $\Psi$ denotes the quark field, and G is the four-fermion interaction coupling constant.

The effective Polyakov-loop potential $\mathcal{U}(\Phi,\bar{\Phi};T)$ that accounts for gauge field self-interaction is a function of the Polyakov-Loop $\Phi$ and its Hermitian conjugate $\bar{\Phi}$

[TABLE]

with the Polyakov line is

[TABLE]

Following Refs. Ruggieri (2011); Gatto and Ruggieri (2012); Cui et al. (2016); Liu et al. (2020), we take the approximation $L^{{\dagger}}=L$ and have $\Phi=\bar{\Phi}$.

The form of the Polyakov effective potential $\mathcal{U}$ is provided by Rößner et al. (2007),

[TABLE]

with

[TABLE]

The parameters are used to reproduce the pure gluonic lattice data, see Table 1. Following Ref. Sakai et al. (2010), we adjust $T_{0}=190\ \mathrm{MeV}$ to take into account the dynamical quarks.

The coupling constant $G$, as pointed out in Refs. Sakai et al. (2010); Kondo (2010), should depend on $\Phi$. Here, we take the form in Ref. Sakai et al. (2010)

[TABLE]

The PNJL model with parameters used in Tables 1 and 2 is in good agreement with lattice QCD data Sakai et al. (2010).

The thermodynamic potential density function $\Omega$ is defined as

[TABLE]

Using finite temperature field theory and mean field approximation, $\Omega$ can be derived as:

[TABLE]

where $M$ is the dynamical quark mass:

[TABLE]

$\sigma=\langle\bar{\Psi}\Psi\rangle$ is the quark chiral condensate. And

[TABLE]

in which $E_{p}=\sqrt{p^{2}+M^{2}}$ is the single quasi-particle energy. The finite temperature contribution term is finite, so here we only impose the cut-off $\Lambda$ on the vacuum term Costa et al. (2010); Ratti et al. (2006); Deb et al. (2009); Fukushima (2004).

Numerical solutions of physical quantities $M$ and $\Phi$ with respect to $\mu$ and $T$ can be obtained by minimizing the thermodynamic potential

[TABLE]

The quantitative study of the QCD phase transition is based on thermal susceptibility

[TABLE]

II.2 Nonextensive PNJL model

The nonextensive version of the thermodynamic potential density function $\Omega_{q}$ is defined as Santos et al. (2014); Ke-Ming Shen and Wang (2017)

[TABLE]

Combining q-algebra and following the same derivation steps as Eq. (II.1), the thermodynamic potential density function $\Omega_{q}$ can be derived as:

[TABLE]

where

[TABLE]

It should be noted that, as a simplification, the form of the Polyakov-loop potential remains unchanged. That is to say, it is only indirectly affected by the nonextensive effects through the saddle point equations. Besides, the parameters of the PNJL model are as usual. Here, we treat $q$ as a thermodynamic variable on the same footing as $T$ and $\mu$. That is to say, we based on the ansatz that the parameters determined at $T=0$, $\mu=0$ and $q=1$ can be used to study the whole region. The same simplifications also appear in Ref. Ke-Ming Shen and Wang (2017).

It is worth noting that to ensure that $e_{q}(x)$ is always a non-negative real function, the following condition must be satisfied:

[TABLE]

During the calculation, we found that the $q$ value cannot be too large, otherwise the CEP position we care about will fall into the non-physical region. Therefore, in this paper we mainly focus on $1\leq q\leq 1.2$. This is also the typical region of $q$ values found by high-energy collisions. Cleymans et al. (2013); Li et al. (2013); Marques et al. (2015); Azmi and Cleymans (2014). In addition, there is no nonextensive effects at low temperature. Because as $T$ tends to zero, $\Omega_{q}$ tends to $\Omega$ as long as $q>1$.

According to Eq. (15), the numerical results for $M$ and $\Phi$ can be obtained by solving the following nonlinear coupling equations:

[TABLE]

where the q-version of the Fermi-Dirac distribution is

[TABLE]

and

[TABLE]

As expected, for $q\rightarrow 1$, the standard distribution function of the usual PNJL model is recovered.

III QCD phase transition and baryon-number susceptibilities within tsallis statistics

III.1 QCD phase transition

The changes of $M$ and $\Phi$ with temperature for different parameters $q$ are shown in Figs. 1, 2. We find that the pseudo-critical temperature decreases with increasing $q$, and the same conclusion can be seen in Ref. Ke-Ming Shen and Wang (2017). Besides, even if the nonextensive effects are taken into account, the pseudo-critical temperature of the chiral and deconfinement crossover transition is still the same. For example, at $q=1$ and $q=1.15$, the pseudo-critical temperature is $173\ \mathrm{MeV}$ and $148\ \mathrm{MeV}$, respectively.

When considering the finite chemical potential, for $q=1$, the variation of the constituent quark mass $M$ with temperature is shown in Fig. 3. It can be clearly seen that the first-order phase transition and the crossover transition occur at chemical potentials $\mu=190\ \mathrm{MeV}$ and $\mu=150\ \mathrm{MeV}$, respectively. The determination of the CEP position is based on the peak of the thermal susceptibility $\chi_{T}$. Because in the CEP position, $\chi_{T}$ is divergent. From Fig. 4, we determine the location of CEP as $(\mu_{c},T_{c})=(171\ \mathrm{MeV},159.5\ \mathrm{MeV})$.

The impact of the nonextensive effects on the location of CEP, which we are most concerned about, is shown in Fig. 5. We find that the nonextensive effects have a significant impact on the CEP position. For $1\leq q\leq 1.08$, with the increase of $q$, the CEP position moves rapidly in an approximate straight line to the direction of lower temperature and larger chemical potential. But for $1.08\leq q\leq 1.2$, we find that the CEP position has an obvious non-monotonic dependence on $q$. Although $T_{c}$ keeps decreasing with the increase of $q$, $\mu_{c}$ first increases and then decreases. The critical value $q_{c}$ is 1.15. The numerical results are shown in Table 3. This non-monotonic behavior also appeared in our previous article Zhao et al. (2021). However, compared with Ref. Zhao et al. (2021), the quantitative difference is significant. This is because for $q=1$, the coupling constant $G(\Phi)$ doubles the value of $T_{c}$, and the nonextensive effects become more pronounced only for sufficiently high temperatures. In addition, another difference is that in Ref. Zhao et al. (2021), as $q$ increase, $T_{c}$ also moves in the opposite direction, and the corresponding critical value $q$ is 1.1. The reappearance of this non-monotonic behavior deserves our attention, and it may be model-independent. The reason why Refs. Rożynek and Wilk (2016); Ke-Ming Shen and Wang (2017) did not find this phenomenon is most likely because they only calculated $q$ to 1.1. Our results are meaningful for search of the CEP position in the relativistic heavy-ion collision, where strong fluctuations and long-range correlation drag the system into out-of-equilibrium states, and BG statistics fails. For a wide class of such systems, it has been shown in recent years that the correct approach is to use Tsallis statistics instead Tirnakli and Borges (2016); Cirto et al. (2018). We find that in the search for CEP, the nonextensive effects cannot be ignored. And the search for CEP in larger or smaller quark chemical potential regions depends on $q_{c}$.

III.2 Baryon-number susceptibilities

In this chapter, we pay attention to the baryon-number susceptibilities because they are related to the moments of the conserved net-baryon number distribution, such as the variance $\sigma^{2}$, the skewness $S$, and the kurtosis $\kappa$. Non-monotonic dependence of these moments, especially high-order quantities, on collision energy $\sqrt{S_{NN}}$ is suggested as an experimental signature of the QCD critical point Stephanov (2009). The detailed formulas are as follows:

[TABLE]

Where the nth order of the baryon-number susceptibility, $\chi_{B}^{(n)}$ is

[TABLE]

According to Eq. (II.2), the quark number density $\rho(T,\mu)$ is

[TABLE]

Here for brevity, we omit the subscript $q$ of the quark number density.

In order to compare with the experiment, one needs the values of $T$ and $\mu_{B}$ along the freeze-out curve as a function of the collision energy $\sqrt{S_{NN}}$. Here we use the empirical fitting formula in Ref. Andronic et al. (2010), as follows

[TABLE]

with the limiting temperature $T_{lim}=164\ \mathrm{MeV}$.

In Fig. 6, $S\sigma$, $\kappa\sigma^{2}$ are shown as a function of $\sqrt{S_{NN}}$. The latest experimental data come from $Au+Au$ collisions at RHIC Adam et al. (2021). By changing the $q$ value ($1\leq q\leq 1.2$, $\Delta q=0.01$), it is not difficult to find that our results are generally consistent with the experimental data at $q=1.07$, especially for $\sqrt{S_{NN}}>19.6\ \mathrm{GeV}$, which fits the experimental data well. But for $\sqrt{S_{NN}}=11.5,14.5\ \mathrm{GeV}$, there is a larger deviation. The results of $\kappa\sigma^{2}$ are qualitatively consistent with Ref. Bzdak et al. (2020). This non-monotonic variation with $\sqrt{S_{NN}}$ shows that the chemical freeze-out points cross the crossover transition line, see Fig. 7. And the intersection of the two is about $\sqrt{S_{NN}}=20\ \mathrm{GeV}$. However, it should be noted that based on our model results, $S\sigma$ also has a non-monotonic dependence on the collision energy, which is qualitatively different from the experimental data and deserves further study. This non-monotonic behavior is also seen in the NJL model Fan et al. (2019). In addition, we also test the dependence of our results on the $T$, $\mu$, taking $T(\sqrt{S_{NN}})\pm 1\ \mathrm{MeV}$ and $\mu(\sqrt{S_{NN}})\pm 4\ \mathrm{MeV}$. From Figs. 6, 7, it can be found that near the intersection, there is a large uncertainty in the results. This is due to a large uncertainty of quark-number susceptibility in this region. All in all, our results generally agree with the experimental data at $q=1.07$, which to some extent indicates that the nonextensive effects are worth considering in relativistic heavy ion collisions. And it may also provide some useful clues for related research to interpret experimental data from a nonextensive perspective.

IV Summary and Conclusion

In this paper, combined with the Tsallis statistics and the PNJL model, first of all, we investigated the sensitivity of the QCD phase transition to deviations from usual BG statistics. At zero chemical potential and finite temperature, we found that the pseudo-critical temperature of the chiral and deconfinement transition decreases as $q$ increases. At finite chemical potential and finite temperature, the most interesting thing is that we found the non-monotonic dependence of CEP on $q$. At the beginning, as $q$ increases, CEP moves rapidly toward a lower temperature and larger quark chemical potential. But when $q\geq 1.15$, CEP turns to move in a direction with lower temperature and quark chemical potential. This means that searching for CEP in larger or smaller quark chemical potential regions in relativistic heavy-ion collisions depends on a critical value $q_{c}$.

Secondly, we studied the moments of the net-baryon number distribution and obtained that $S\sigma$, $\kappa\sigma^{2}$ are generally consistent with the experimental data when $q=1.07$. This indicates that the nonextensive effects may be worth considering in relativistic heavy-ion collisions. In addition, quark stars, as candidates for observed massive stars ($\geq 2M_{\odot}$), have attracted much attention in astronomy Li et al. (2020); Chu et al. (2019); Li et al. (2019); Chen et al. (2016). Therefore, studying the influence of nonextensive effects on the structure and evolution of protoquark stars is a very meaningful topic Lavagno and Pigato (2011).

Finally, how to introduce the nonextensive effects directly into the pure gluonic sector is a problem worthy of further study. As mentioned above, in the standard Polyakov-loop potential, there is no space to introduce nonextensive effects because this potential is based on group integral rather than momentum integral. Can we try to consider other gluonic potentials, for example the Meisinger-Miller-Ogilvie model Meisinger et al. (2002), which can describe the confinement-deconfinement nature of QCD, qualitatively. Of course, this will be a brand new attempt, and these issues are our future research directions.

Acknowledgements.

This work is supported by the National Natural Science Foundation of China (under Grants No. 12005192), and the Project funded by China Postdoctoral Science Foundation (Grant No. 2020M672255 and No. 2020TQ0287).