Critical Exponents of the Chiral Phase Transition in the PNJL Model

It is a significant topic to explore the phase transition of strongly interacting matter in both theory and heavy-ion collision experiment. The calculations of lattice QCD^[1–6] indicate that, at high temperature and small chemical potential, the transformation from quark-gluon plasma to hadrons is a smooth crossover. Although the simulation at large chemical potential is still not available in lattice QCD due to the sign problem, some effective quark models which incorporate the fundamental symmetries of QCD (e.g., Refs. [7–19]) indicate the existence of a first-order phase transition at large chemical potential, with a critical end point (CEP) connecting with the crossover transition.

Searching for the CEP is just one of the primary tasks of RHIC^[20-21], and some possible indictions related to the critical phenomenon were discovered in the first phase of beam energy scan (BES I)^[22–24]. To pin down the existence of CEP, the second phase of the beam energy scan and the fixed target experiments will be performed at RHIC STAR. Besides, other planned experiments to probe the phase transition at lower temperatures and larger chemical potentials will be done at FAIR/NICA/J-PARC. The progress in experiments will provide us good opportunities to identify the QCD phase structure.

Experimentally, the event-by-event measurements of fluctuations of conserved charges undertake the mission of searching for the QCD phase structure^{[20-21, 25–36]}. In the beam energy scan experiments, the high-order cumulants of net proton (proxy for net baryon) have been measured, and the multiplicity distributions, skewness and kurtosis, are exploited to identify the QCD phase transition. When approaching to the vicinity of CEP, the density fluctuations and the corresponding multiplicity distribution grow fast and diverge at the CEP^[37]. The divergent strength is correlated to the critical exponents. Therefore, the study of the critical exponents is meaningful for the analysis of particle multiplicity distributions and identifying the location of CEP.

Although some critical behaviors of chiral transition have been investigated in literature (e.g., Refs. [38–44]), the critical exponents related to multiplicity distributions of conserved charges such as net-baryon number skewness $ S\sigma $ and kurtosis $ \kappa\sigma^2 $ are still not investigated, which are more closely connected to the experimental measurements. Calculating the critical exponents of baryon number skewness and kurtosis is one of the main tasks in this work.

On the other hand, for the chiral transition at low density and large chemical potential, the idealized first-order phase transition is usually derived based on the equilibrium thermodynamics in literature (e.g., Refs. [7–19, 37, 45-46]). The associated nonequilibrium effects of created matter is ignored, however, for an expanding system, the nonequilibrium effect possibly plays a crucial role on the phase transformation. In heavy-ion collision experiments, the unstable and/or metastable phase of chiral transition may be triggered with the collision energy is reduced to a few GeV^{[44, 47–55]}. At the spinodal boundary where the unstable phase and metastable phase are separated, strongly density fluctuations will generate. The corresponding multiplicity distribution of conserved charges will also diverge, which possibly bring new disturbance to identify the first-order phase transition.

Considering the planned experiments with reduced collision energies to explore the critical phenomenon and the first-order transition as well as the observables taken in measurement, we will, in this study, investigate the critical exponents of baryon number skewness and kurtosis at both the CEP and the spinodal boundary. This investigation is meaningful for the dada analysis to identify the critical behavior and the first-order chiral transition.

The Lagrangian density in the three flavor PNJL (Poyakov-Nambu–Jona-Lasinio) model is taken as

where $ q $ denotes the quark fields with three flavors, $ u,\ d $, and $ s $; $ \hat{m}_{0} = {\texttt{diag}}(m_{u},\ m_{d},\ m_{s}) $ in flavor space; $ G $ and $ K $ are the four-point and six-point interacting constants, respectively. The $ \hat{\mu} = diag(\mu_u,\mu_d,\mu_s) $ are the quark chemical potentials.

The covariant derivative in the Lagrangian is defined as $ D_\mu = \partial _\mu-iA_\mu $. The gluon background field $ A_\mu = \delta_\mu^0A_0 $ is supposed to be homogeneous and static, with $ A_0 = g{\cal{A}}_0^\alpha \frac{\lambda^\alpha}{2} $, where $ \frac{\lambda^\alpha}{2} $ is $ SU(3) $ color generators. The effective potential $ {\cal{U}}(\Phi[A],\,\bar{\Phi}[A],\,T) $ is expressed in terms of the traced Polyakov loop $ \Phi = ({\rm{Tr}}_c L)/N_C $ and its conjugate $ \bar{\Phi} = ({\rm{Tr}}_c L^\dagger)/N_C $. The Polyakov loop $ L $ is a matrix in color space

where $ \beta = 1/T $ is the inverse of temperature and $ A_4 = {\rm{i}}A_0 $.

The Polyakov-loop effective potential used in this study is

where $ a(T) = a_0+a_1(\frac{T_0}{T})+a_2(\frac{T_0}{T})^2 $ and $ b(T) = b_3(\frac{T_0}{T})^3 $. The parameters $ a_0 = 3.51 $, $ a_1 = -2.47 $, $ a_3 = 15.2 $, and $ b_3 = -1.75 $ were derived in Ref. [56] by fitting the thermodynamics of pure gauge sector in LQCD. $ T_0 = 210 $ MeV is implemented when fermion fields are included.

In the mean field approximation, the constituent quark masses $ M_i $ are obtained as

where $ \phi_i $ stands for quark condensate of the flavor $ i $. The thermodynamical potential is derived as

where $ E_i = \sqrt{p^{2}+M_i^2} $ is the energy-momentum dispersion relation.

The values of $ \phi_u, \phi_d, \phi_s, \Phi $ and $ \bar{\Phi} $ are determined by minimizing the thermodynamical potential

In the calculation a cut-off $ \Lambda $ is implemented in 3-momentum space for divergent integrations. The model parameters given in Ref. [57] are taken with $ \Lambda = 603.2 $ MeV, $ G\Lambda^{2} = 1.835 $, $ K\Lambda^{5} = 12.36 $, $ m_{u,d} = 5.5 $ and $ m_{s} = 140.7 $ MeV, determined by fitting $ f_{\pi} = 92.4 $ MeV, $ M_{\pi} = 135.0 $ MeV, $ m_{K} = 497.7 $ MeV and $ m_{\eta} = 957.8 $ MeV.

In the PNJL model with the parameters given above, the critical end point locates $ \mu_c = 291$ MeV and $ T_{c} = 131$ MeV. First of all, we investigate the four standard thermodynamical critical exponents $ \alpha,\, \beta, \,\gamma $ and $ \delta $ usually used to describe the critical behavior. Fig. 1 shows the paths to calculate the four critical exponents on the phase diagram.

Figure 1.  (color online) First-order phase transition line with the CEP and the paths taken to calculate the critical exponents.

The exponent $ \alpha $ is defined by the power scaling of specific heat $ C_{\rho} = - T\left(\frac{\partial ^2 \Omega}{\partial T^2}\right)_{\rho} $ along the first-order axis towards the critical endpoint

For the convenience of numerical calculation, the calculation is usually performed along the extension line of the first-order transition. In the numerical analysis, we define $ t = \frac{T-T_{c}}{T_{c}} $, which measure the distance from the critical temperature $ T_{c} $. We plot $ \ln(C_{\rho} $) as a function of $ \ln(t) $ in Fig. 2. From the figure we can see that all points fall on the line which is read $ y = 0.02x+2.11 $ with the best fitting. The slope of the line gives the critical exponent $ \alpha = 0.02 $

Figure 2.  (color online) Logarithm of specific heat approaching to the critical end point along the extension of first-order transition line.

The exponent $ \beta $ describes the discontinuity of quark number density $ \rho $ as the temperature approaches the critical temperature $ T_c $ along the first-order transition line

$ \Delta\rho $ is the difference of quark density on the two sides of the first-order transition for a given temperature. The best fitting of $ \ln(\Delta \rho) $ as a function of $ \ln(t) $ given in Fig. 3 indicates the critical exponent $ \beta = 0.50 $.

Figure 3.  (color online) Logarithm of $ \Delta\rho$ as function of $ \ln(t)$ along the first-order line approaching to the CEP.

The exponent $ \gamma $ is related to the baryon number susceptibility $ \chi_2\equiv\left(\frac{\partial \rho}{\partial \mu}\right)_{_T} $ along the first-order axis

We still take the logarithm on both sides of the above equation. The numerical result of the best fitting is plotted in Fig. 4, which gives the critical exponent $ \gamma = 1.03 $.

Figure 4.  (color online) Logarithm the susceptibility $ \chi_{2}$ approaching to the CEP along extension line of the first-order transition.

The exponent $ \delta $ is defined by the the relationship between $ \rho-\rho_c $ and $ \mu-\mu_c $ when the temperature is fixed to $ T_c $

With the same method used above the logarithm of $ \rho-\rho_{c} $ as a function $ \ln(\mu-\mu_{c}) $ is plotted in Fig. 5. The best fitting shows the critical exponent $ \delta = 3.03 $.

Figure 5.  (color online) Logarithm of $\rho-\rho_c$ and $\mu-\mu_{c}$ as $\mu$ approaches $\mu_{c}$ at the critical temperature.

The numerical analysis indicates that the values of $ \alpha, \beta, \gamma $ and $ \delta $ for strongly interacting matter in the PNJL model also approximately satisfy the following relations:

and

They are consistent with the conclusion in Refs. [58-59] and the analysis of the Landau-Ginzburg theory in the mean-field approximation with $ \alpha = 0, \,\beta = 1/2, \,\gamma = 1 $ and $ \delta = 3 $(Refs. [43, 60] and reference therein).

Besides the four critical exponents above, we also calculate the divergence of $ \chi_2^B $, $ S\sigma $ and $ \kappa \sigma^2 $ at the critical temperature $ T_c $ with $ \mu $ approaching to $ \mu_c $ from the high density side. These critical exponents are defined as

and

respectively. We should note that here $ \nu $ and $ \eta $ are not the critical exponents relevant to correlation length. $ S\sigma $ and $ \kappa \sigma^2 $ are the ratios of high-order susceptibilities: $ S\sigma = \chi^{B}_{3}/\chi^{B}_{2 } $ and $ \kappa \sigma^2 = \chi^{B}_{4}/\chi^{B}_{2} $. The baryon number susceptibility can be derived with

For more details, one can refer to Ref. [37].

The numerical values of these critical exponents are $ \nu = 0.70, \zeta = 1.03 $ and $ \eta = 2.08 $. The two critical exponents $ \zeta = 1.03 $ and $ \eta = 2.08 $ are closely related to baryon number fluctuation distributions, skewness and kurtosis. The collision energy dependence of skewness and kurtosis are taken to analyze the phase transition signatures in the beam energy scan experiments of heavy-ion collision, especially to pin down the critical endpoint. We can see that the critical exponents $ \zeta $ is smaller than $ \eta $. It means that the kurtosis is more sensitive than skewness in the vicinity of CEP. The values of $ \nu $, $ \zeta $ and $ \eta $ will be compared with the corresponding divergence along the spinodal boundary in the next subsection, which are possibly useful to analyze the first-order transition.

In this subsection, we discuss the noneuqilibrium effect of the first-order transition below the critical temperature. Generally, a first-order transition is always associated with the nonequilibrium effect during the condensation in terms of the convex structure of thermodynamic pressure. To a certain degree, the nonequilibrium effect can manifest in a rapid expanding system. Some signals different from the idealized first-order transition will probably be observed in the upcoming HIC(Heavy-ion Collision) experiments at RHIC/FARI/NICA. Therefore, it is meaningful to investigate the density fluctuations and the corresponding multiplicity distributions resulting from the nonequilibrium first-order transition with the inclusion of the unstable (spinodal) and metastable phases. Here we focus on the critical behavior of multiplicity distribution of baryon number fluctuation along the spinodal boundary.

For the convenience of later discussion, we plot the full $ T-\rho^{}_B $ phase diagram of the first-order transition in Fig. 6. The lines marked $ A $ and $ D $ are the binodal boundaries of the idealized first-order transition. This figure demonstrates that the unstable (spinodal) region filled with gray color is surrounded by the metastable regions (two yellow bands). With the increase of temperature, the unstable and metastable regions gradually decrease. Finally, the spinodal and binodal curves meet at the critical end point.

Figure 6.  (color online) Phase structure of the first-order transition including the metastable and unstable regions. $B$ and $C$ are the spinodal boundaies. $A$ and $D$ are the binodal curves of the first-order transition.

As an example, we fix the temperature $ T = 100$ MeV, and show the relations of $ \mu_q, \chi_2^B $, skewness and kurtosis as functions of baryon density in Fig. 7. It is clear that $ \chi_2^B $, skewness and kurtosis diverge at the spinodal boundaries $ B $ and $ C $. We will calculate the critical exponents at the point $ B $ from the low density side and at the point $ C $ from the high-density side.

Figure 7.  (color online) $\mu_q$, $\chi^{B}_{2}$, skewness and kurtosis as functions of baryon density at $T = 100$ MeV.

Near the singularity the behavior of $ \chi_2^B $, skewness and kurtosis can be described by the power law forms

and

where $ \mu_{0} $ is the chemical potential on spinodal boundary ($ B $ or $ C $). In the calculation, we still use the linear logarithmic fitting, thus the Eq. (17)$ \thicksim (19)$ can be written as

and

The critical exponent $ \zeta^\prime $ of skewness along the spinodal lines plotted in Fig. 6 at several temperatures in the first-order transition region are given in Fig. 8. The numerical results shows that $ \zeta^\prime \sim 1.02\pm 0.02 $ for all temperatures^① . The critical exponent $ \eta^\prime $ of kurtosis at different temperatures are given in Fig. 9. The numerical results shows that $ \zeta^\prime \sim 2.08\pm0.02 $ for all temperatures. Comparing with critical exponents $ \eta $ and $ \zeta $ at the CEP, we find that $ S\sigma $ ($ \kappa\sigma^2 $) on the spinodal boundary has the same divergent strength as that at the CEP.

Figure 8.  (color online) Logarithm of skewness with respect to $\ln(|\mu-\mu_{c}|)$ along the isothermal spinodal lines at $T = 80,100,120$, and $131\,$MeV ($T_{c}$). The numerical value is the slope of the line, i.e., the value of $\zeta^\prime$.

Figure 9.  (color online) Logarithm of kurtosis with respect to $\ln(|\mu-\mu_{c}|)$ along the isothermal spinodal lines at $T = 80,100,120$, and $131\,$MeV ($T_{c}$). The numerical value is the slope of the line, i.e., the value of $\eta^\prime$.

For the critical exponent $ \nu' $, it increases gradually along the spinodal line from lower temperatures to higher ones. The evolving behaviors of $ \nu' $ along the spinodal boundary are plotted in Fig. 10. The numerical calculation show that $ \nu' = 0.52\pm0.02 $ for $ T<110 $ MeV. A rapid increase occurs near the critical region, and $ \nu'\approx0.7 $ (0.66) when approaching the CEP along line $ C $ ($ B $) in Fig. 6. These results are consistent with the two-flavor NJL model, which give $ \nu' = 1/2 $ at spinodals and $ \nu' = 2/3 $ at the CEP^[50]. In the chiral limit, $ \nu'\sim 1/2 $ at both the spinodal line and CEP, but for finite quark mass, the analysis in the Ginsburg-Landau theory shows that $ \nu'\sim 2/3 $ at the CEP and 1/2 at the spinodal line far from the CEP. The smooth evolution of the singularity from the spinodal lines to the CEP is expected and is confirmed by the numerical result of the PNJL model. The evolution of $ \nu' $ toward the CEP also indicates a change of universality class.

Figure 10.  (color online) Values of $\nu'$ at different temperatures along the spinodal boundaries. The red and blue lines are the results along the line $B$ and $C$ in Fig. 6, respectively.

In the following, we discuss some possible phase transition signals in observation. For the idealized first-order chiral phase transition in equilibrium, the divergence of density fluctuations only takes place at the critical endpoint and then predicts a peak structure of kurtosis as a function of collision energy. However, when the metastable phase is considered, the calculations show that the divergence of density fluctuations can also occur at the spinodal boundary of the first-order phase transition. Therefore, some unique signals of the chiral phase transition, different from those idealized phase transition in equilibrium, may be generated in experiments. For baryon number kurtosis, since the critical exponent along the spinodal line is the same as that at the CEP, the predicted peak structure as a function of collision energy in equilibrium thermodynamics may be indistinct, and kurtosis always takes relatively larger values in the first-order phase transition region. It weakens the distinction between the critical fluctuations and fluctuations from the first-order phase transition.

In this research, we investigated the critical exponents of chiral transition in the PNJL model. The four standard critical exponents, $ \alpha, \beta, \gamma $ and $ \delta $, are consistent with Landau-Ginzburg theory in the mean-field approximation. Besides, we calculated the criticality of skewness and kurtosis of net-baryon number fluctuations. The critical exponent $ \eta $ of kurtosis is larger than the critical exponent $ \zeta $ of skewness at the CEP. It means that the measurement of kurtosis is more sensitive than skewness if the critical region can be reached in heavy-ion collision.

Considering the spinodal instability may be trigged in heavy-ion collisions with the collision energies decreased to a few GeV, we also studied the critical behavior of skewness, kurtosis and $ \chi_2^B $ along the spinodal boundaries. We found that the critical exponent of skewness (kurtosis) has the same divergent strength as that at the CEP, but the critical exponent $ \nu^\prime $ has a rapid increase when approaching to the CEP along the spinodal boundary. Since the divergence of skewness and kurtosis at the spinodal boundaries, the signals to identify the first-order transition in the future experiments will be disturbed to a certain degree. Some deviations from the prediction of standard first-order transition may be found in measurement.