Introduction

The search for a first-order phase transition and critical point (CP) in the Quantum Chromodynamics (QCD) phase diagram is a fundamental goal of relativistic heavy-ion collision experiments [1-4]. Lattice QCD (LQCD) simulations have shown that the phase transition between the quark-gluon plasma (QGP) and hadronic matter is a smooth crossover at nearly zero baryon chemical potential (μ_B～0) [5-7]. Based on investigations from various effective models, such as the Nambu-Jona-Lasinio (NJL) model and quark-meson (QM) model, as well as advanced functional methods including the Dyson-Schwinger equation (DSE) and functional renormalization group (FRG) [8-20], the transition can be of the first-order in the large μ_B region, resulting in a critical end point on the first-order phase transition line. However, different theoretical models and methods exhibit significant differences in the predicted CP location, even for the same model with different parameter sets [21]. To further explore the QCD phase diagram and search for a possible CP signal, experimental programs such as the Beam-Energy Scan (BES) are currently underway at the Relativistic Heavy-Ion Collider (RHIC). These involve varying the collision energy of heavy-ion collisions to cover a wide range of temperatures T and μ_B in the QCD phase diagram [22-24]. In the BES program, physical observations sensitive to CP and/or first-order phase transitions have been made, such as the pion Hanbury-Brown Twiss (HBT) radii [25, 26], baryon directional flow [27, 28], net proton fluctuation [29, 30], light nuclear yield ratio [31-33], and interaction of charged hadrons [34]. Although nonmonotonic dependencies of these observations on collision energy have been observed, significant uncertainties remain in the nonmonotonic energy range.

Fluctuations in the conserved charges such as net charge, net baryon number, and net strangeness are predicted to depend on the nonequilibrium correlation length (ξ), which diverges at the singular critical point in the ideal system. However, they are typically limited by the finite size or the finite time of the system, and thus serve as possible signals of CP of the strongly interacting QCD matter created in relativistic heavy-ion collisions [35, 36]. Recent reports from first phase of the beam energy scan program (BES-I) of the RHIC [29, 30] demonstrated the non-monotonic variation of high-order moments of proton multiplicity (replacing the net baryon number) with collision energy in the range $7.7\text{ GeV}<\sqrt{s_{\text{NN}}}<27\text{ GeV}$. The measured higher-order moments refer to the skewness ($S=\langle {\left(\delta N\right)}^3\rangle /{\sigma }^3$) and kurtosis ($\kappa =\langle {\left(\delta N\right)}^4\rangle /{\sigma }^4-3$), where $\delta N=N-M$, M is the mean and σ is the quadratic variance. Furthermore, the cumulant ratio, $C_4/C_2=\kappa {\sigma }^2$, of proton multiplicity distribution from $\sqrt{s_{\text{NN}}}=3.0\text{ GeV}$ Au + Au collisions was reported [24]. Compared with the Poisson baseline, a significant suppression in protons C₄/C₂=-0.85 was observed in the most central 0%-5% collisions at 3 GeV. This can be explained by fluctuations driven by baryon number conservation in a region of high baryon density, where hadronic interactions are dominant. These results imply that the QCD critical region, if created in heavy-ion collisions, can only exist at energies higher than 3 GeV. Light nuclei production has been suggested for exploring the existence and location of the QCD first-order phase-transition region with a corresponding critical point. In the 0%-10% most central Au + Au collisions at $\sqrt{s_{\text{NN}}}=19.6\text{ GeV}$ and 27 GeV, $N_{\text{t}}\times N_{\text{p}}/N_{\text{d}}^2$ shows enhancements relative to the coalescence baseline, with a combined significance of 4.1σ [33]. Further studies on heavy-ion collisions are required to confirm whether these nonmonotonic dependencies are due to large fluctuations near the CP and to narrow the energy range at which CP may occur.

Hydrodynamics simulations provide another option for exploring QCD phase transitions. The spatiotemporal evolution of QCD matter can be successfully described using relativistic dissipative hydrodynamics [37-39]. The speed of sound is a crucial physical quantity in hydrodynamics, which carries important information for describing the evolution of strongly interacting matter and final observables in heavy-ion collision experiments. According to Ref. [40-42], the speed of sound as a function of the charged particle multiplicity $\langle \text{d}{N}_{\text{ch}}/\text{d}\eta \rangle$ can be extracted from heavy-ion collision data. A recent study Ref. [43] estimated the value of c_s as well as its logarithmic derivative with respect to the baryon number density and attempted to build a connection with the cumulants of the net baryon number in matter created in heavy-ion collisions to aid in detecting the QCD CP. Furthermore, the behavior of cs as a function of the baryon number density influences the equation of state (EOS) of neutron star matter. An EOS with a pronounced peak in the speed of sound was found to leave a clear and unique signature in the main frequency of the postmerger gravitational wave (GW) spectrum [44], which can sensitively probe the hadron-quark phase transition in the dense core. However, the polytropic index, defined as $\gamma \equiv \partial \ln P/\partial \ln \epsilon$, is potentially important for determining the presence of quark matter in the core of massive neutron stars [45-48]. The approximate rule that γ < 1.75 [45] or γ ≤ 1.6 and $c_{\text{s}}^2\le 0.7$ [49] can be used as a good criterion for separating quark matter from hadronic matter in a massive neutron star core. A recent study [50] introduced a new conformality criterion quantity $d_{\text{c}}\equiv \sqrt{{\text{Δ}}^2+{\left({\text{Δ}}^{\text{'}}\right)}^2}$, where $\text{Δ}\equiv (\epsilon -3P)/(3\epsilon )$ and ${\text{Δ}}^{\text{'}}\equiv \text{d}\text{Δ}/\text{d}$ ln ε are the normalized trace anomaly and its logarithmic rate of change with respect to the energy density, respectively, and compared d_c with the speed of sound, polytropic index, and normalized trace anomaly. The results indicate that the conformality criterion d_c<0.2 is considerably more restrictive than the criterion γ <1.75 for quark matter in compact stars. Therefore, the speed of sound and polytropic index can serve as new probes for QCD CP or first-order phase transitions in future experimental exploration.

We utilized the significant findings of phase diagrams and critical points in the experimental and theoretical research mentioned above in recent years to explore the properties of QCD phase diagrams in the present work by calculating high-order moments, speed of sound, and polytropic index in an extended 3-flavor NJL model including an eight-quark scalar-vector interaction. Scalar-vector interactions are important for reproducing nuclear saturation properties when using the NJL-type model for nuclear matter [51]. When applied to the hadron-quark phase transition, the strength of the scalar-vector interactions can change the critical temperature T_c, thereby facilitating the study of the effect of CP in heavy-ion collisions [52, 53]. The remainder of this paper is organized as follows. In Sect. 2, we briefly review the calculation process of the extended 3-flavor NJL model, as well as the calculation formulas for the high-order moments, speed of sound, and polytropic index. The phase diagram of the high-order moments, speed of sound, and polytropic index, as well as the effects of the scalar-vector interaction on the phase diagram, are discussed in Sect. 3. Finally, the summary is presented in Sect. 4.

Theoretical model

The Lagrangian density of the extended three-flavor NJL model with an eight-quark scalar-vector interaction is given by [53] $\begin{array}{c}{\mathcal{L}}_{\text{NJL}}^{SU\left(3\right)}=\overline{\psi }\left(i{\gamma }^{\mu }{\partial }_{\mu }-\widehat{m}\right)\psi +G_{\text{S}}{\displaystyle \sum _{a=0}^8\left[{\left(\overline{\psi }{\lambda }^a\psi \right)}^2+{\left(\overline{\psi }i{\gamma }_5{\lambda }^a\psi \right)}^2\right]}\\ -K\left[\det \overline{\psi }\left(1+{\gamma }_5\right)\psi +\det \overline{\psi }\left(1-{\gamma }_5\right)\psi \right]\\ +G_{\text{SV}}｛{\displaystyle \sum _{a=1}^3\left[{\left(\overline{\psi }{\lambda }^a\psi \right)}^2+{\left(\overline{\psi }i{\gamma }_5{\lambda }^a\psi \right)}^2\right]}｝\\ \times ｛{\displaystyle \sum _{a=1}^3\left[{\left(\overline{\psi }{\lambda }^{\mu }{\lambda }^a\psi \right)}^2+{\left(\overline{\psi }{\gamma }_5{\lambda }^{\mu }{\lambda }^a\psi \right)}^2\right]}｝\mathrm{,}\end{array}$ (1) where $\psi ={\left(u\mathrm{,}d\mathrm{,}s\right)}^{\text{T}}$ represents the three flavor quark fields and $\widehat{m}=diag\left(m_{\text{u}}\mathrm{,}{m}_{\text{d}}\mathrm{,}{m}_{\text{s}}\right)$ is the corresponding current quark mass matrix. In the above formula, ${\lambda }^a\left(a=\mathrm{1,}\dots \mathrm{,8}\right)$, where ${\lambda }^0=\sqrt{2/3}$ are the Gell-Mann matrices; G_S and G_SV are the strengths of the scalar coupling and scalar-vector coupling, respectively; and the K term is the six-quark interaction (also known as the Kobayashi-Maskawa-t’Hooft interaction) that breaks the axial U_A(1) symmetry [54, 55]. This term gives rise to six-point interactions in the three flavors and is responsible for the flavor mixing effect.

In the mean-field approximation [56], quarks can be treated as quasiparticles whose constituent mass Mi is determined by the gap equation $\begin{array}{c}{M}_{\text{u}}=m_{\text{u}}-2G_{\text{S}}{\varphi }_{\text{u}}+2K{\varphi }_{\text{d}}{\varphi }_{\text{s}}\\ -2G_{\text{SV}}{\left({\rho }_{\text{u}}+{\rho }_{\text{d}}\right)}^2\left({\varphi }_{\text{u}}+{\varphi }_{\text{d}}\right)\mathrm{,}\\ M_{\text{d}}=m_{\text{d}}-2G_{\text{S}}{\varphi }_{\text{d}}+2K{\varphi }_{\text{u}}{\varphi }_{\text{s}}\\ -2G_{\text{SV}}{\left({\rho }_{\text{u}}+{\rho }_{\text{d}}\right)}^2\left({\varphi }_{\text{u}}+{\varphi }_{\text{d}}\right)\mathrm{,}\\ M_{\text{s}}=m_{\text{s}}-2G_{\text{S}}{\varphi }_{\text{s}}+2K{\varphi }_{\text{u}}{\varphi }_{\text{d}}\mathrm{,}\end{array}$ (2) where ${\varphi }_i=\langle \overline{i}i\rangle$ and ρi denote the quark condensate and net quark number densities, respectively. The thermodynamic properties of the three-flavor quark matter are determined by the partition function $\mathcal{Z}=Tr\left[\exp \left[-\beta \left(\widehat{H}-\mu \widehat{N}\right)\right]\right]$. Here, β =1/T and μ are the reciprocals of the temperature and chemical potential, respectively, while $\widehat{H}$ and $\widehat{N}$ represent the Hamiltonian and the corresponding conserved charge number operators, respectively. The thermodynamic potential of the system can be calculated as $\begin{array}{c}{\text{Ω}}_{\text{NJL}}^{SU(3)}=G_{\text{S}}\left({\varphi }_{\text{u}}^2+{\varphi }_{\text{d}}^2+{\varphi }_{\text{s}}^2\right)+3G_{\text{SV}}{\left({\varphi }_{\text{u}}+{\varphi }_{\text{d}}\right)}^2{\left({\rho }_{\text{u}}+{\rho }_{\text{d}}\right)}^2\\ -4K{\varphi }_{\text{u}}{\varphi }_{\text{d}}{\varphi }_{\text{s}}-2N_{\text{c}}{\displaystyle \sum _{i=u\mathrm{,}d\mathrm{,}s}{\displaystyle {\int }_0^{\Lambda }}}\frac{{\text{d}}^{3}p}{{\left(2\pi \right)}^3}{E}_i\\ -2T{\displaystyle \sum _{i=u\mathrm{,}d\mathrm{,}s}{\displaystyle {\int }_0^{\Lambda }}}\frac{{\text{d}}^{3}p}{{\left(2\pi \right)}^3}\left(z^{+}\left(E_i\right)+z^{-}\left(E_i\right)\right)\mathrm{,}\end{array}$ (3) where N_c =3 is the number of colors and $z^{\pm }\left(E_i\right)=N_{\text{c}}\ln \left[1+e^{-\beta \left(E_i\mp {\tilde{\mu }}_i\right)}\right]\mathrm{,}$ (4) with an effective chemical potential $\begin{array}{l}{\tilde{\mu }}_{\text{u}}={\mu }_{\text{u}}+2G_{\text{SV}}\left({\rho }_{\text{u}}+{\rho }_{\text{d}}\right){\left({\varphi }_{\text{u}}+{\varphi }_{\text{d}}\right)}^2\mathrm{,}\\ {\tilde{\mu }}_{\text{d}}={\mu }_{\text{d}}+2G_{\text{SV}}\left({\rho }_{\text{u}}+{\rho }_{\text{d}}\right){\left({\varphi }_{\text{u}}+{\varphi }_{\text{d}}\right)}^2\mathrm{,}\\ {\tilde{\mu }}_{\text{s}}={\mu }_{\text{s}}\mathrm{.}\end{array}$ (5) The energy Ei is expressed as $E_i=\sqrt{M_i^2+p^2}$, where Mi is the mass of the constituent quarks. The quark condensate ${\varphi }_i$ and net quark number density ρi can be derived by minimizing the thermodynamic potential as follows: $\begin{array}{l}{\varphi }_i=2N_{\text{c}}{\displaystyle \int }\frac{{\text{d}}^{3}p}{{\left(2\pi \right)}^3}\frac{M_i}{E_i}\left(f_i+{\overline{f}}_i-1\right)\mathrm{,}\\ {\rho }_i=2N_{\text{c}}{\displaystyle \int }\frac{{\text{d}}^{3}p}{{\left(2\pi \right)}^3}\left(f_i-{\overline{f}}_i\right)\mathrm{,}\end{array}$ (6) where $\begin{array}{l}{f}_i=\frac{1}{1+e^{\beta \left(E_i-{\tilde{\mu }}_i\right)}}\mathrm{,}\\ {\overline{f}}_i=\frac{1}{1+e^{\beta \left(E_i+{\tilde{\mu }}_i\right)}}\mathrm{.}\end{array}$ (7) In the present numerical calculation, we employ the following values for the parameters [57] as $m_{\text{u}}=m_{\text{d}}=5.5\text{ MeV}$, $m_{\text{s}}=135.7\text{ MeV}$, $G_{\text{S}}{\Lambda }^2=3.67$, $K{\Lambda }^5=9.29$, and the cutoff value in the momentum integral $\Lambda =631.4\text{ MeV}$. For the coupling constant G_SV, Eqs. (2) and (5) show that its value can affect the effective masses of the quarks and their effective chemical potentials. However, the entire term, including G_SV does not affect the description of QCD vacuum properties at zero baryon density. Furthermore, treating the value of G_SV as a free parameter allows us to obtain the different properties of quark matter.

In the following section, we consider the fluctuation moments of conserved quantities, such as the net baryon number, from the above 3-flavor NJL model. The experimentally measured proton multiplicity (net baryon number) distributions are described by the moments or ratios of the cumulants. The relationships between the cumulants Cn and moments are defined as $C_1=M$, $C_2={\sigma }^2$, $C_3=S{\sigma }^3$, and $C_4=\kappa {\sigma }^4$. Here, M and σ are the mean and quadratic variance, respectively, while the skewness S and kurtosis κ refer to higher-order moments. The ratios of the cumulants are often used to reduce the volume dependence, which can be read as $C_2/C_1={\sigma }^2/M$, $C_3/C_2=S\sigma$, and $C_4/C_2=\kappa {\sigma }^2$. Moreover, the ratios of cumulants in the theoretical calculations are related to the ratios of the baryon-number susceptibilities. The nth-order susceptibility of the net baryon number in the grand ensemble can be expressed as the derivative of the thermodynamic potential [58] ${\chi }_n^{\text{B}}=\frac{{\partial }^n\left(-\text{Ω}/T\right)}{\partial {\left({\mu }_{\text{B}}/T\right)}^n}\mathrm{.}$ (8) The ratios of the baryon-number susceptibilities related to the skewness Sσ and kurtosis κσ² measured experimentally in relativistic heavy-ion collisions can be calculated using [58] $S\sigma \left(B\right)=\frac{C_3^{\text{B}}}{C_2^{\text{B}}}=\frac{{\chi }_3^{\text{B}}}{{\chi }_2^{\text{B}}}\mathrm{,}\kappa {\sigma }^2\left(B\right)=\frac{C_4^{\text{B}}}{C_2^{\text{B}}}=\frac{{\chi }_4^{\text{B}}}{{\chi }_2^{\text{B}}}\mathrm{.}$ (9)

In heavy-ion collision experiments, the speed of sound c_s is a crucial thermodynamic observable, which provides important information for describing the evolution of the fireball. Recently, the authors in Ref. [43] attempted to build a connection between c_s and the cumulants of the baryon number distribution in matter created in heavy-ion collisions to aid in detecting QCD CP. The general definition of the speed of sound is $c_{\text{s}}^2=\partial P/\partial \epsilon$, and its calculation requires specifying the thermodynamic variables that are kept constant. For the QGP created in relativistic heavy-ion collisions, it evolves with a constant entropy per baryon s/ρ_B, which justifies using the adiabatic speed of sound $c_{s/{\rho }_{\text{B}}}$. In contrast, the isothermal speed of sound, c_T, is widely used in neutron star matter. However, the values of $c_{s/{\rho }_{\text{B}}}^2$ and $c_{\text{T}}^2$ largely coincide near the critical point [59, 60], implying that the cumulants of the baryon number distribution in heavy-ion collisions can be used to estimate the isothermal squared speed of sound [43]. Therefore, in this study, we mainly discuss the phase diagram properties of the isothermal speed of sound. The pressure and entropy density in the NJL model can be derived using the thermodynamic relations in the grand canonical ensemble as follows: $P=-{\text{Ω}}_{\text{NJL}}\mathrm{,}s=-\frac{\partial {\text{Ω}}_{\text{NJL}}}{\partial T}\mathrm{,}$ (10) and the energy density can be calculated as $\epsilon =-P+Ts+{\mu }_{\text{B}}{\rho }_{\text{B}}\mathrm{.}$ (11) Using the Jacobi method and thermodynamic relations, $c_{\text{T}}^2$ can be written in terms of T and μ_B as $c_{\text{T}}^2={\left(\frac{\partial P}{\partial \epsilon }\right)}_{\text{T}}=\frac{{\rho }_{\text{B}}}{T{\left(\frac{\partial s}{\partial {\mu }_{\text{B}}}\right)}_{\text{T}}+{\mu }_{\text{B}}{\left(\frac{\partial {\rho }_{\text{B}}}{\partial {\mu }_{\text{B}}}\right)}_{\text{T}}}\mathrm{.}$ (12) Meanwhile, the polytropic index introduced in [45] $\gamma \equiv \partial \ln P/\partial \ln \epsilon$ serves as a criterion for separating hadronic and quark matter [45, 47-49]. The isothermal polytropic index can be derived using the expression for the speed of sound: ${\gamma }_{\text{T}}={\left(\frac{\partial P}{\partial \epsilon }\right)}_{\text{T}}/\frac{P}{\epsilon }=\frac{\epsilon }{P}{c}_{\text{T}}^2\mathrm{.}$ (13) The above expression shows that the speed of sound is an important intermediate quantity for calculating the polytropic index. The differences between these two quantities may lead to a polytropic index that includes more intricate features near the phase boundary, which could potentially provide a new probe for exploring QCD CP.

Results and Discussions

First, we discuss the higher-order susceptibilities of the net baryon number in the μ_B-T plane based on the 3-flavor NJL model with different scalar-vector coupling constants G_SV. As shown in Fig. 1, we present the skewness Sσ(B) and kurtosis κσ²(B) of the baryon number fluctuations as functions of temperature and baryon chemical potential. The red and blue regions represent positive and negative values, respectively, and the green regions represent the values of Sσ(B) and κσ²(B) that are approximately zero. The black dashed-dotted and solid lines represent the chiral crossover and first-order phase transition, respectively, while the white dots connecting the chiral crossover and first-order phase transition represent the CPs. The chiral phase transition boundaries separate the red and blue regions for skewness, whereas the phase boundaries for kurtosis pass through the blue regions and divide the red regions with the CP located at the ends of the blue regions. In addition, the phase transition boundaries and CPs vary depending on the scalar-vector coupling constants G_SV. Although decreasing the value of G_SV effectively enhances the temperature of the CP, its baryon chemical potential decreases. The effects of G_SV can be understood using Eq. (2) and (5). According to Eq. (5), the negative G_SV resembles a vector interaction in the NJL model, which induces a repulsive interaction among the quarks. Meanwhile, compared with the scalar term G_S in the NJL model, a negative G_SV counteracts the repulsive interaction, as shown in Eq. (2), which reduces the constituent quark mass. For quark matter at a low baryon density (or chemical potential), the repulsive quark interaction with a negative G_SV is stronger than the attractive quark interaction [53, 57], resulting in a net repulsive effect of the negative G_SV. For quark matter at intermediate densities, because the quadratic of G_SV term in Eq. (2) depends on the quark density, the effect of the G_SV term on constituent quark masses is significantly enhanced with increasing quark density, which generates an effectively attractive interaction among quarks. Consequently, when the coupling constant $G_{\text{SV}}=-200{\Lambda }^{-8}$, the CP temperature increases to 104 MeV, while the chemical potential decreases to 750 MeV. For a positive G_SV, the effects on the properties of quark matter are opposite to those of a negative G_SV leading to a lower temperature and higher baryon chemical potential CP. Therefore, the range of the CP location obtained by varying G_SV is sufficiently large to cover the region that can be probed during heavy-ion collisions. In quark matter at very high densities, where the chiral symmetry is largely restored and the quark condensates are close to zero, the effects of the G_SV term become less important. This differs from the usual vector interaction in the NJL model that becomes stronger at high densities.

Fig. 1Full-screen View

(Color online) Skewness Sσ(B) and kurtosis κσ²(B) of the baryon number fluctuations as functions of the temperature and baryon chemical potential in the NJL model with different scalar-vector coupling constants G_SV. The black dash-dotted and solid lines represent the chiral crossover and first-order phase transition, respectively, while the white dots connecting the chiral crossover and first-order phase transition represent the CPs

In relativistic heavy-ion collision experiments, fluctuations in the net baryon number and net charge number were measured at the chemical freeze-out lines. However, the location of chemical freeze-out could not be accurately determined at the RHIC-BES energies. Several criteria exist for chemical freeze-out, such as fixed energy per particle at approximately 1 GeV, fixed total density of baryons and anti-baryons, fixed entropy density over T³, and the percolation model. Recently, the hypothetical freeze-out lines derived by fitting the experimental data have been expressed as $T=a-b{\mu }_{\text{B}}^2-c{\mu }_{\text{B}}^4$ [2, 61]. However, the exact correspondence between the experimental parameters and location on the phase diagram remains unclear. Therefore, for simplicity and better elucidation of the properties of the phase diagram near the critical point, the hypothetical freeze-out lines shown in Fig. 2 were obtained by rescaling μ_B of the phase boundaries with factors of 0.90, 0.95, 0.98, 1.02, and 1.05, which can also be found in Refs. [16, 20, 62]. Figure 3 shows the skewness Sσ(B) and kurtosis κσ²(B) of the net baryon number fluctuations as functions of the baryon chemical potential along the hypothetical chemical freeze-out lines with different scalar-vector coupling constants G_SV. Along the hypothetical chemical freeze-out lines, nonmonotonic dependencies in both Sσ(B) and κσ²(B) occur, and peak or dip structures appear around the critical points. In contrast, κσ²(B) features both a positive peak and negative dip, whereas Sσ(B) has either a single peak or dip, whose nature (peak or dip) depends on whether the hypothetical freeze-out lines are below the chiral phase transition boundaries. The amplitude of the peak or dip structures increases significantly as the hypothetical chemical freeze-out lines approach the CPs, for example, along the hypothetical lines with factors of 0.98 and 1.02. The nonmonotonic behavior of κσ²(B) is qualitatively consistent with the experimental results, while the experimental results for Sσ(B) seem to be positive above $\sqrt{s_{\text{NN}}}=7.7\text{ GeV}$ from both STAR and PHENIX measurements [27]. Further experiments on collisions of different heavy ions at different energies are required, so that the freeze-out lines are close to or even above the phase transition boundaries. In addition, the peak and dip structures of the net-baryon susceptibilities shift to the low-baryon chemical potential (and high-temperature) side with a negative scalar-vector coupling constant G_SV, which is consistent with our previous discussions and expectations. Furthermore, the amplitudes of variation for both Sσ(B) and κσ²(B) values significantly increase as G_SV decreases, which may be because G_SV is a quartic term in the thermodynamic potential. Thus, fluctuations in the density and quark condensation near the CP are greatly amplified.

Fig. 2Full-screen View

(Color online) Hypothetical chemical freeze-out lines with different scalar-vector coupling constants G_SV by rescaling μ_B of the chiral phase transition boundary with factors of 0.90, 0.95, 0.98, 1.02, and 1.05

Fig. 3Full-screen View

(Color online) Skewness Sσ(B) and kurtosis κσ²(B) of the net baryon number fluctuations as functions of the baryon chemical potential along the hypothetical chemical freeze-out lines with different scalar-vector coupling constants G_SV

The speed of sound, which is a fundamental property of a substance, provides important information for describing the evolution of strongly interacting matter and the final observables in heavy-ion collision experiments. It can be obtained from heavy-ion collision data through the charged particle multiplicity $\langle d{N}_{\text{ch}}/d\eta \rangle$ [40-42]. Furthermore, the net baryon fluctuations in heavy-ion collisions can be used to estimate the isothermal speed of sound and its logarithmic derivative [43]. The numerical results from effective models [59, 60, 67] suggest that the speed of sound at the CP is the global minimum of the full phase diagram and is also considered a sensitive probe of the hadron-quark phase transition in the compact star core. Fig. 4 presents the contour maps of the squared speed of sound $c_{\text{T}}^2$ in the μ_B-T plane based on the NJL model with different scalar-vector coupling constants G_SV. The values of $c_{\text{T}}^2$ in the chiral breaking regions corresponding to the blue regions are mostly less than 0.2. For the different scalar-vector coupling constants G_SV, changes in the range of chiral breaking (blue) regions and the location of the critical points are evident. After passing through the chiral phase transition boundary, $c_{\text{T}}^2$ rapidly increases with the restoration of chiral symmetry. In Fig. 4, a dip structure appears on the low-temperature side of the chiral phase transition boundary, and a minimum of $c_{\text{T}}^2$ appears around the critical point. For the different scalar-vector coupling constants, the value of $c_{\text{T}}^2$ is less than 1 in the entire μ_B-T plane, satisfying the constraint of causality. However, in the NJL model, $c_{\text{T}}^2$ does not approach the value $c_{\text{s}}^2=1/3$ of conformal matter limit at the high temperatures and baryon chemical potentials, because of the cutoff parameter set for the momentum integral in the NJL model. For the Polyakov-looped Nambu-Jona-Lasinio (PNJL) model, because momentum integration is finite and does not require an ultraviolet cutoff [11, 12], the squared speed of sound at high temperatures and baryon chemical potentials approaches 1/3 [60]. To further explore the properties of the speed of sound near the critical points, we present the results of the squared speed of sound $c_{\text{T}}^2$ along the hypothetical chemical freeze-out lines with different scalar-vector coupling constants G_SV in Fig. 5. For the hypothetical freeze-out lines over the chiral phase transition boundaries, the $c_{\text{T}}^2$ curves rapidly increase after softening near the CPs. However, the speed of sound along the hypothetical freeze-out lines below the chiral phase transition boundary varies considerably: $c_{\text{T}}^2$ rapidly decreases near the CPs, especially as it approaches the local minimum along the curve by a factor of 0.98, followed by a spinodal behavior that eventually continues to decrease. Additionally, the spinodal regions are influenced by G_SV. As G_SV decreases, the span of the spinodal regions expands accordingly; conversely, it decreases or approaches zero with increasing G_SV. This is also observed in Fig. 1, where the positive G_SV reduces the amplitude of the fluctuations around the critical points.

Fig. 4Full-screen View

(Color online) Contour maps of the squared speed of sound $c_{\text{T}}^2$ in the μ_B-T plane based on the NJL model with different scalar-vector coupling constants G_SV. The black dash-dotted and solid lines indicate the chiral crossover and first-order phase transition, respectively, while the white dots connecting the chiral crossover and first-order phase transition represent the critical points (CPs)

Fig. 5Full-screen View

(Color online) Squared speed of sound $c_{\text{T}}^2$ as a function of baryon chemical potential along the hypothetical chemical freeze-out lines with different scalar-vector coupling constants G_SV

Subsequently, we investigated the contour maps of the polytropic index γ_T in the μ_B-T plane with different G_SV values, as shown in Fig. 6. The values of γ_T are almost less than 1.5 in the chiral breaking regions with different G_SV. On either side of the phase boundary, particularly around the first-order phase transition, non-monotonic behavior, including dip (blue region) and peak (red region) structures in γ_T are clearly observed. Following an isotherm that passes through a first-order phase transition, the polytropic index rapidly decreases to a minimum value, increases rapidly to the maximum with chiral restoration, and subsequently decreases again, eventually stabilizing at high temperatures and baryon chemical potentials. This can be explained by Eq. (13), where the polytropic index is equal to the squared speed of sound divided by $P/\epsilon$. Around the first-order phase transition, the squared speed of sound (the derivative of P with respect to ε) and $P/\epsilon$ are not synchronized; therefore, a local minimum and maximum of the polytropic index appear at different chemical potentials. A more detailed explanation is provided in Ref. [60]. In Fig. 6, the amplitudes and spans of the dip and peak structures are influenced by the coupling constants G_SV. The amplitude and span of the peak structures (red regions) are effectively enhanced by decreasing the value of G_SV; however, those of the dip structures (blue regions) decrease. Similar to the speed of sound, we investigate the properties of γ_T along the hypothetical chemical freeze-out lines in Fig. 7. All curves of γT first decrease and attain a minimum near the CP, and closer critical points correspond to smaller minimum values, such as the hypothetical freeze-out lines with 0.98 or 1.02. Consistent with the results for Sσ(B) and κσ²(B), the mininum value of γ_T shifts to a low baryon chemical potential with a negative G_SV. Further, a negative G_SV leads to a larger minimum value of γ_T because the negative G_SV, which serves as a repulsive interaction at a low chemical potential, stiffens the equation of stateand increases the speed of sound and polytropic index. Compared to $c_{\text{T}}^2$ as shown in Fig. 5, the polytropic index exhibits a more pronounced dip structure near the CP, which could provide a more sensitive probe for the QCD CP in future experiments.

Fig. 6Full-screen View

(Color online) Contour maps of the polytropic index γ_T in the μ_B-T plane based on the NJL model with different scalar-vector coupling constants G_SV. The black dash-dotted and solid represent the chiral crossover and first-order phase transition, respectively, and the white dots connecting the chiral crossover and first-order phase transition represent the critical points (CPs)

Fig. 7Full-screen View

(Color online) Polytropic index γ_T as a function of the baryon chemical potential along the hypothetical chemical freeze-out lines for different scalar-vector coupling constants G_SV

Summary and Outlook

In conclusion, we investigated the properties of the phase diagrams of higher-order susceptibilities, speed of sound, and polytropic index based on an extended NJL model with an eight-quark scalar-vector interaction. We found that all the aforementioned quantities exhibited non-monotonic behavior around the phase transition boundaries, which also revealed the location of the critical point. Specifically, the baryon susceptibility κσ²(B) featured both a positive peak and negative dip, whereas Sσ(B) had either a single peak or dip, whose nature depended on whether the hypothetical freeze-out lines were below the phase-transition boundaries. Meanwhile, for $c_{\text{T}}^2$, a dip structure appeared on the low-temperature side of the phase transition boundary, and a minimum appeared around the critical point. Along the hypothetical freeze-out lines below the phase transition boundary, the speed of sound rapidly decreased near the CP, followed by spinodal behavior, whereas the polytropic index γ_T exhibited a more pronounced dip as it approached the CP. Compared to the speed of sound, the polytropic index could provide a more sensitive probe for QCD CP in future experiments. In this study, we explored the effect of the eight-quark scalar-vector interaction on the phase diagram. The chiral phase transition boundary and critical point varied depending on the scalar-vector coupling constants G_SV. At finite densities and temperatures, the negative G_SV term exhibited attractive interactions, which enhanced the critical point temperature and reduced the chemical potential. The G_SV term also affected the high-order susceptibilities, speed of sound, and polytropic index. The nonmonotonic (peak or dip) structures of these quantities shifted to a low baryon chemical potential (and high temperature) with a negative G_SV. G_SV also changed the amplitude and range of the nonmonotonic regions. Therefore, the scalar-vector interaction is useful for locating the phase boundary and critical point in the QCD phase diagram by comparing the experimental data. Although significant uncertainty remains regarding the properties of the QCD phase diagram and the precise location of the critical point, the study of the nonmonotonic behavior of these observables is of great interest. A more accurate BES-II measurement in the near future will provide further information regarding the QCD phase diagram. Moreover, in the study of massive neutron stars, the approximate rule γ < 1.75 [45] or γ ≤ 1.6 and $c_{\text{s}}^2\le 0.7$ may be used as a good criterion for separating quark matter from hadronic matter in a massive neutron star core [49]. Recently, a new conformality criterion d_c<0.2, defined as comprising γ and c_s, was introduced and compared with the speed of sound, polytropic index, and normalized trace anomaly [50]. The results indicate that the conformality criterion d_c<0.2 is considerably more restrictive than the criterion γ <1.75 for quark matter in compact stars [50]. In recent years, machine learning (ML)and high-energy nuclear physics have begun to merge, yielding interesting results [67]. ML methods have also been applied to phase transition studies, enabling finer characterization of the QCD phase diagram and deeper understanding of quark matter properties [68]. Therefore, it is worth looking forward to new methods and observations related to high-order susceptibilities, speed of sound, and polytropic index and their application to the search for critical points in heavy-ion collisions and the study of compact stars.