Nonlinear susceptibilities of heavy ion collisions within the Polyakov-loop-extended Nambu-Jona-Lasinio model

In recent decades, the quantum chromodynamics (QCD) phase transition and corresponding phase diagram have been studied extensively through a series of theoretical tools and heavy ion collision experiments [1]. In particular, QCD transitions from hadronic matter to the quark-gluon plasma (QGP) is considered occur in the heavy ion collision experiments [2]. This process is a smooth crossover at small baryon chemical potential and high temperature, and a first-order phase transition is expected at high baryon chemical region. The intersection between first-order phase transition and crossover region is called QCD critical end point (CEP), people want to find CEP by experiment.

From the usual view, quark-number (or baryon-number) susceptibility (the second order) will show some singularity [3, 4] near the CEP. Many phenomenological models [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17] have been used and lattice QCD [18, 19, 20] calculations are used to find the location of the CEP. Baryon number fluctuations have long been considered to be closely related to phase transitions, experiments such as STAR’s beam energy scan (BES) program and PHENIX experiments at the Relativistic Heavy-Ion Collider (RHIC) [21] have been carried out to measure them. In our work, we used the experimental data of RHIC, which obtained the fluctuation of baryon number in Au + Au collisions. The $n$th cumulant of baryon-number fluctuations is proportional to the $n$th order of baryon-number susceptibilities [22, 2, 23], so we need to calculate the nth order of baryon-number susceptibilities theoretically in order to compare them with the experimental values.

It can be found that the quark-number density determined by the corresponding dressed quark propagator at finite chemical potential [24]. There are variational methods that one can use a bare propagator and obtain nonperturbative results, such as hard thermal loop perturbation theory, optimized perturbation theory, and so no. However, we don’t use variational methods. Here, we will get the dressed quark propagator through the PNJL model [25], which is the Nambu-Jona-Lasinio (NJL) model improved with the Polyakov loop. In recent years, the PNJL model has been successfully used to describe the thermodynamics of QCD with two and two-plus-one flavors, it allows for a simultaneous computation of quantities sensible to confinement and chiral symmetry breaking [26]. Some other models such as NJL model, the chiral random matrix model, the linear sigma model and chiral perturbation theory are based on chiral symmetry, but they all lack any dynamics coming form the Polyakov loop. [27].

The rest of this paper is organized as follows: In Sec.II nonlinear susceptibilities and a mean field description of the PNJL model is presented. In Sec.III we present our calculation results, which were transformed into experimentally observable results for comparisons and research. Finally, in Sec.IV we will summarize our results and give the conclusions.

The renormalized partition function of QCD at zero temperature and finite chemical is of the form [24]

	$\displaystyle\mathcal{Z}[\mu]=$	$\displaystyle\int\mathcal{D}\overline{q}_{R}\mathcal{D}q_{R}\mathcal{D}A_{R}\exp\{-S_{R}\left[\overline{q}_{R},q_{R},A_{R}\right]+$		(1)
		$\displaystyle\int d^{4}x\mu Z_{2}\overline{q}_{R}(x)\gamma_{4}q_{R}(x)\}$

where $S_{R}\left[\overline{q}_{R},q_{R},A_{R}\right]$ is standard renormalized Euclidean QCD action, $q_{R}$ is the renormalized quark field with three flavors and three colors; $Z_{2}=Z_{2}\left(\zeta^{2},\Lambda^{2}\right)$ is quark wave-function renormalization constant, where $\zeta$ is the renormalization point and $\Lambda$ is the regularization mass-scale. By leave the ghost field term and its integration measure to be understood, the pressure density $\mathcal{P}(\mu)$ is given by

$$\mathcal{P}(\mu)=\frac{1}{\mathcal{V}}\ln\mathcal{Z}[\mu]$$

(2)

where $\frac{1}{\mathcal{V}}$ is the four-volume normalizing factor. And the quark-number density is

$$\rho(\mu)=\frac{\partial\mathcal{P}(\mu)}{\partial\mu}$$

(3)

From Eq.(1), (2) and (3), we can easily obtain the result [28]

$$\rho(\mu)=(-)N_{c}N_{f}Z_{2}\int\frac{d^{4}p}{(2\pi)^{4}}\operatorname{tr}\left[G[\mu]({p})\gamma_{4}\right]$$

(4)

Where $N_{c}$ and $N_{f}$ are, respectively, the number of colors and flavors; $G[\mu]({p})$ is the quark propagator; From this formula, we can see that in the case of zero temperature and limited chemical potential, the quark number density is determined only by the dressed quark propagator at finite chemical potential, by ignore the $\mu$ dependence of the dressed gluon propagator and assume that the dressed quark propagator at finite $\mu$ is analytic in the neighborhood of $\mu$=0, then we can obtain the following expression [29]

$$G[\mu]({p})^{-1}=G(\tilde{p})^{-1}$$

(5)

Where $G[\mu]({p})^{-1}=i\gamma$ · ${p}+M$, $M$ is effective quark mass; $\tilde{p}=\left(\vec{p},p_{4}+i\mu\right)$. By some mathematical methods, like the matsubara frequency, we can get the expression for the quark-number density at finite temperature

$$\rho(T,\mu)=(-)N_{c}N_{f}T\sum_{i=-\infty}^{+\infty}\int\frac{d^{3}\vec{p}}{(2\pi)^{3}}\operatorname{tr}\left[G\left(\tilde{p}_{n}\right)\gamma_{4}\right]$$

(6)

the fourth component of momentum is $\omega_{n}+i\mu$, and the fermion frequencies $\omega_{n}=(2n+1)\pi T$. A baryon consists of three quarks, so the quark-number density is three times as many as the baryon-number density. Then the $\rho_{B}$ to the baryon chemical potential $\mu_{b}$’s (n-1) order derivatives are defined by the nonlinear susceptibilities of baryons of order n [30].

	$\displaystyle\chi_{B}^{(n)}$	$\displaystyle=\frac{\partial^{n-1}}{\partial\mu_{B}^{n-1}}\rho_{B}=\frac{\partial^{n-1}}{3^{n}\partial\mu^{n-1}}\rho(T,\mu)$		(7)
		$\displaystyle=(-)\frac{N_{c}N_{f}T}{3^{n}}\sum_{i=-\infty}^{+\infty}\int\frac{d^{3}\vec{p}}{(2\pi)^{3}}\operatorname{tr}_{\gamma}\left[\frac{\partial^{n-1}G\left(\tilde{p}_{n}\right)}{\partial\mu^{n-1}}\gamma_{4}\right]$

And then by taking the trace and summing the frequencies, we can get

		$\displaystyle\rho(T,\mu)=$		(8)
		$\displaystyle(-)\frac{N_{c}N_{f}T}{\pi}\int(\frac{1}{e^{\frac{\sqrt{M^{2}+\vec{p}^{2}}-\mu}{T}}+1}-\frac{1}{e^{\frac{\sqrt{M^{2}+\vec{p}^{2}}+\mu}{T}}+1})p^{2}dp$

So what we’re going to do is figure out the effective quark mass $M$ in the PNJL model.

The Lagrangian density of two flavors of equal-mass quarks in the PNJL model is

	$\displaystyle\mathcal{L}=$	$\displaystyle\overline{q}(\gamma\cdot D+m)q-G\left[(\overline{q}q)^{2}+\left(\overline{q}i\gamma_{5}\tau q\right)^{2}\right]$		(9)
		$\displaystyle+\mathcal{U}(\Phi,\overline{\Phi};T)$

where $m$ is the common current-quark mass; $D_{\mu}={\partial_{\mu}}+iA_{\mu}$, with $A_{\mu}(x)=g_{s}A_{\mu}^{a}\lambda^{a}/2$ describe the matrix-valued gluon field configuration that fits the model; G is the four-fermion interaction strength and $\mathcal{U}$ is a Polyakov-loop effective potential.

Adopting the mean-field approximation, we take $L=\overline{L}$ from the beginning as in [31], the effective potential of the model can be separately expressed as [32, 33],

$$\Phi=\frac{1}{N_{c}}\operatorname{Tr}_{c}L=\overline{\Phi}$$

(10)

according to the classical background field in Eq.(11), the effective potential of Polyakov rings is expressed as follows [34]

	$\displaystyle\frac{1}{T^{4}}\mathcal{U}(\overline{\Phi},\Phi;T)=$	$\displaystyle\frac{1}{T^{4}}\mathcal{U}(\Phi;T)$		(11)
	$\displaystyle=$	$\displaystyle-\frac{1}{2}a(T)\Phi^{2}+b(T)$
		$\displaystyle\times\ln\left[1-6\Phi^{2}+8\Phi^{3}-3\Phi^{4}\right]$

with $(\overline{t}=T_{0}/T)$

$$a(\overline{t})=a_{0}+a_{1}\overline{t}+a_{2}\overline{t}^{2},\quad b(\overline{t})=b_{3}\overline{t}$$

(12)

For the sake of reproduce the lattice results for pure-gauge QCD chromodynamics and the T-dependence of Polyakov loop, we set these parameters and list them in Table I [34].

The four-fermion coupling, $G$, is considered a constant [27], The effective potential of the PNJL model is obtained by means of the mean-field approximation [35, 36],

	$\displaystyle\Omega=$	$\displaystyle\Omega\left(M,\Phi;T,\mu\right)$		(13)
	$\displaystyle=$	$\displaystyle\mathcal{U}(\Phi;T)+\frac{(M-m)^{2}}{4G}-4N_{c}\int\frac{d^{3}\vec{p}}{(2\pi)^{3}}\omega$
		$\displaystyle-8T\int\frac{d^{3}\vec{p}}{(2\pi)^{3}}\ln\left[\mathcal{F}\right]$

where $M$ and $\mathcal{F}$ is

$$\omega=\sqrt{M^{2}+\vec{p}^{2}}$$

(14)

$$\mathcal{F}_{\pm}=1+3\Phi\left[\mathrm{e}^{-\frac{\omega_{\pm}}{T}}+\mathrm{e}^{-2\frac{\omega_{\pm}}{T}}\right]+\mathrm{e}^{-3\frac{\omega_{\pm}}{T}}$$

(15)

$\omega_{\pm}=\omega\pm\mu$. It’s worth noting that, we have not yet explicitly solved the regularization problem of the PNJL model. The last term in the second line of Eq. (13) is a divergent quantity, a regularization procedure must be introduced. For the reasons explained in [33], we impose a hard cutoff on both the integrals, although, the last one is not divergent. We employ a hard cutoff $\Lambda$, which is listed in table I . Using that value, and $m$ and $g$ listed, we obtained a good description of in-vacuum pion properties. At this point, we can determine the quark mass gap evolution with intensive parameters by solving two external conditions simultaneously:

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

(16)

We solved the above equation by iterative method, then obtained the high-order susceptibility by numerical differentiation method.

In order to compare our results with the experimental date and other methods, we mainly calculate the following formulas [22]

	$\displaystyle S\sigma$	$\displaystyle=\frac{T\chi_{B}^{(3)}}{\chi_{B}^{(2)}}$		(17)
	$\displaystyle\kappa\sigma^{2}$	$\displaystyle=\frac{T^{2}\chi_{B}^{(4)}}{\chi_{B}^{(2)}}$
	$\displaystyle\frac{\kappa\sigma}{S}$	$\displaystyle=\frac{T\chi_{B}^{(4)}}{\chi_{B}^{(3)}}$

where $\sigma^{2}$ is the variance, $S$ is the skewness, and $\kappa$ is the kurtosis. In order to better compare with other researchers’ work, we show $S\sigma$, $\kappa\sigma^{2}$ and $\frac{\kappa\sigma}{S}$ as a function of $\sqrt{S_{NN}}$ for $Au+Au$ collisions at RHIC in Fig.1. The top column represents the corresponding freeze-out chemical potential $\mu_{B}$ to $\sqrt{S_{NN}}$. Here we use the more commonly used empirical relationship to show the correlations between $\sqrt{S_{NN}}$ and the bulk properties ($\mu_{B}$ and $T$) of chemical freeze-out [37, 38, 39, 40].

	$\displaystyle T\left(\mu_{B}\right)$	$\displaystyle=a-b\mu_{B}^{2}-c\mu_{B}^{4}$		(18)
	$\displaystyle\mu_{B}(\sqrt{S_{NN}})$	$\displaystyle=\frac{d}{1+e\sqrt{S_{NN}}}$
	$\displaystyle\mu_{B}$	$\displaystyle=3\mu$

where $a=0.166\pm 0.002$ GeV, $b=0.139\pm 0.016$ $\text{GeV}^{-1}$, $c=0.053\pm 0.021$ $\text{GeV}^{-3}$, $d=1.308\pm 0.028$ GeV, and $e=0.273\pm 0.008$ $\text{GeV}^{-1}$ [39]. According to the above formula, we list the corresponding values of $T,\mu_{B},\mu$, and $\sqrt{S_{NN}}$ in Table II.

The value of $\sqrt{S_{NN}}$ is for comparison with experimental data.The results that we obtained by the PNJL model are compared with lattice QCD [41], DSEs [42] and experimental data [21].

As can be seen from Fig.1, the fitting of experiment of the PNJL model in $S\sigma$ and $\kappa\sigma^{2}$ is not as good as that of lattice QCD. However, we have a better fit at $\frac{\kappa\sigma}{S}$. This result is very similar to that of DSEs, so we used a similar approach to explore the reasons for the inconsistency with the experiment. We fixed the chemical potential $\mu$ at 9, 27 and 76 $MeV$, and we compute the curves of $S\sigma$, $\kappa\sigma^{2}$, and $\frac{\kappa\sigma}{S}$ as a function of temperature $T$, from 100 to 160 $MeV$.

We present the curves of $S\sigma$, $\kappa\sigma^{2}$, and $\frac{\kappa\sigma}{S}$ obtained by our calculation as functions of $T$ and $\mu$ in Fig.2 , where $\mu$ = 9, 27 and 76 $MeV$ respectively. We found that $S\sigma$ and $\kappa\sigma^{2}$ in our results were smaller than the experimental values, and $\frac{\kappa\sigma}{S}$ is larger than the experimental value. DSEs also encountered the same problem, fix $T$ at 160 and 166 MeV, $S\sigma$ and $\kappa\sigma^{2}$ obtain by DSEs are all too small to compare with the experimental data when $\mu$ is less than 90 MeV. Therefore, we did the same calculation by the PNJL model, the results are shown in Fig.3. We find that there is an intersection between $S\sigma$ and the minimum experimental value, but it is smaller than other experimental values. And $\kappa\sigma^{2}$ is always smaller than experimental value no matter how $\mu$ changes.

Lowering the value of $m$, or $\Lambda$, or increasing the value of $G$ will bring our calculation closer to the experimental value. However, the parameters of the model are all obtained by fitting the meson mass $m_{\pi}$, the pion decay constant $f_{\pi}$ and the quark condensates $\langle\overline{q}q\rangle$ [43], so we don’t know how they should change with temperature and chemical potential. In this paper, we didn’t fine-tune the parameters.

In this paper, we calculate baryon-number susceptibilities from the Polyakov-loop-extended Nambu-Jona-Lasinio models, then compare our results with the results from the heavy ion collision experiments and other models.

The curves of $S\sigma$, $\kappa\sigma^{2}$, and $\frac{\kappa\sigma}{S}$ obtained by our calculation as functions of $T$ and $\mu$ are very different from the DSEs curves in areas outside the scope of the experiment, but the results of the two different models were surprisingly consistent at several points in the experimental values. In addition, $S\sigma$ and $\kappa\sigma^{2}$ of DSEs are smaller than the experimental value no matter how $\mu$ changes, and our results are similar to those of DSEs. This means that the different in DSEs and the PNJL model will lead to different peak point and CEP calculated by the two methods, and that is also why the curves of $S\sigma$, $\kappa\sigma^{2}$, and $\frac{\kappa\sigma}{S}$ as function of $T$ and $\mu$ obtained by the two methods are different. The consistency of the two different methods at the experimental point shows the validity and the rationality of the two methods in dealing with the heavy ion collision experiments. However, $S\sigma$ and $\kappa\sigma^{2}$ are smaller than the experimental results, which means that the PNJL model also has the same characteristics as DSEs, and there are still further reasons for us to explore. In general, the PNJL model is more concise and convenient than DSEs, but it still gives reasonable results. This proves that the PNJL model can be well used to deal with this issue. The difference between the results in Fig.2 and DSE also provides the correct direction for the experiment to verify the two models. Outside the existing experimental area, the two methods give significantly different results, so further research is needed.

Obviously, there’s more that can be done here, such as our Polyakov loop potential is $\mu$-independent and we don’t consider the impact of a vector interaction, $G_{V}(\overline{q}\gamma_{\mu}q)^{2}$, which may influence the results at high baryonic densities. There are some other work where Polyakov loop potenial is $\mu$-dependent [44] or $G_{V}(\overline{q}\gamma_{\mu}q)^{2}$ [45] is considered, our subsequent work will dive into these issues.

The author would like to thank Dong Bai for helpful communications. This work is supported by the National Natural Science Foundation of China (Grant No. 11535004, 11761161001, 11375086, 11120101005, 11175085, 11235001 and 11975167).