Exploring the equation of state (EOS) and the phase structure of quark gluon plasma (QGP) is a hot topic in nuclear physics^[1−15]. The speed of sound in QGP is a crucial mechanical characteristic to describe the variation of EOS. Its dependence on temperature and baryon density carries important information in describing the evolution of a fireball and the final observables. Recently, there have been some attempts to extract the speed of sound of QGP from the experimental data. For example, the studies indicate that the speed of sound is a function of charged particle multiplicity $ \langle {\rm{d}}N_{{\rm{ch}}}/{\rm{d}}\eta \rangle $^[16−17]. An interesting method to estimate the speed of sound is proposed by building its connection with the net baryon number cumulants and the quantum chromodynamics (QCD) phase structure^[18].

Besides, there has been growing interest in the speed of sound in neutron star matter^[19−21]. The variation of the speed of sound with baryon density affects the mass-radius relation and the tidal deformability of neutron stars, and provides a useful probe for studying the EOS of neutron star matter. To obtain a neutron star with a mass greater than twice that of the sun, some studies show that the EOS of neutron star matter must be sufficiently stiff in some density range, corresponding to a speed of sound squared much greater than 1/3^[22−28]. The speed of sound in neutron star matter is also connected to the gravitational wave frequencies induced by the g-mode oscillation of a neutron star^[29]. Furthermore, the magnitude of speed of sound affects the dynamics of primordial density perturbations in the era of cosmic QCD phase transition, which is relevant for the induced gravitational wave signals^[30].

Some calculations about the speed of sound in QCD matter have been performed in lattice QCD^[31−35], (P)NJL model^[36−40], quark-meson coupling model^{[7, 41]}, hadron resonance gas (HRG) model^[42−43], field correlator method (FCM)^[44−45] and quasiparticle model^[46]. In particular, the relationship between the speed of sound and the QCD phase transition has been intensively studied^[47−48]. The numerical results indicate that the dependence of the speed of sound on the system parameters (temperature, density and chemical potential) is closely related to the phase structure of QCD.

By far, the studies in the literature mainly focus on the speed of sound in symmetric quark matter, and there is a lack of systematic research on the speed of sound in the isospin asymmetric quark matter. In this work we will explore the behavior of the speed of sound in asymmetric quark matter. On the one hand, we will derive the formulae for calculating the speed of sound in asymmetric quark matter with two conserved charges. On the other hand, we will demonstrate numerically the dependence of the speed of sound on isospin asymmetric parameter in the PNJL model. The relevant research has certain significance in exploring the equation of state of QCD and the evolution of QGP.

The general definition of speed of sound is

where $ p $ and $ \epsilon $ are pressure and energy density. A specifying constant quantity X is required to describe the propagation of the compression wave through a medium. For a fireball created in relativistic heavy-ion collisions, it evolves with a constant entropy density to baryon number density ratio $s/\rho_{\rm{B}}^{}$ under the ideal fluid approximation. Therefore, the speed of sound can be calculated along the isentropic curve

In practice the above definition can only be used to calculate the speed of sound along isentropic trajectories. To calculate the speed of sound throughout the full phase diagram, it is necessary to derive the corresponding formula as a function of temperature T and baryon number density $\rho_{\rm B}^{}$. With the fundamental thermodynamic relation, Eq. (2) can be rewritten as

where $\sigma=s / \rho_{{\rm{B}}}^{}$ is defined for convenience.

The speed of sound will be calculated under the physical condition that the number density of strange quark is zero $ \, \rho_{\rm s}=0 $, which is the requirement of strangeness conservation in strong interaction. The asymmetric parameter $ \alpha $ is defined as the ratio of u, d quark number density $ \alpha=\rho_{\rm{u}}/\rho_{\rm{d}} $. Then, we can derive the following formulae in the grand canonical ensemble

In the above formulae, $ \mu_{\rm u} $, $\mu_{\rm d}^{}$ represent the chemical potential of u, d quark. For a given asymmetric parameter $ \alpha $, Eq. (3) can be used to calculate the speed of sound in the $T-\rho_{\rm B}^{}$ plane with these formulae. For the case of $ \alpha=1 $, these formulae return to those for symmetric quark matter.

In the numerical calculation, we take the 2+1 flavor PNJL quark model. the Lagrangian density is written as

where $ q $ denotes the quark fields with three flavors, $ u,\ d $, and s; $\hat{m}_{0}={\rm{diag}}(m_{{\rm{u}}},\; m_{{\rm{d}}},\; m_{{\rm{s}}})$ in flavor space; $ G $ and $ K $ are the four-point and six-point interacting constants, respectively. The $\hat{\mu}={\rm{diag}}(\mu_{\rm {\rm{u}}},\; \mu_{\rm d},\;\mu_{\rm{s}})$ are the quark chemical potentials. The covariant derivative in the Lagrangian is defined as $D_\mu=\partial_\mu - {\rm{i}}A_\mu$. The gluon background field $ A_\mu=\delta_\mu^0A_0 $ is supposed to be homogeneous and static, with $A_0=g{\mathcal{A}}_0^\alpha \frac{\lambda^\alpha}{2}$, where $ \frac{\lambda^\alpha}{2} $ is the ${\rm S} {\rm U}(3)$ color generators. The effective potential $U(\varPhi[A],\;\bar{\varPhi}[A],\;T)$ is expressed with the traced Polyakov loop $\varPhi=({\rm{Tr}}_{\rm{c}} L)/N_{\rm{C}}$ and its conjugate $\bar{\varPhi}=({\rm{Tr}}_{\rm{c}} L^\dagger )/N_{\rm{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 is

where

The parameters $ a_0=3.51 $, $ a_1=-2.47 $, $ a_2=15.2 $, $ b_3=-1.75 $ and $ T_0=270 $ MeV are fitted according to the earlier lattice simulation of QCD thermodynamics in the pure gauge sector^[49]. When the quark fields are included $ T_0=210 $ MeV will be implemented to describe well the thermodynamics and phase transition of QCD. The constituent quark mass in the mean field approximation can be derived as

where $ \phi_i $ stands for quark condensate of the flavor $ i $.

The thermodynamic potential of bulk quark matter is derived as

where ${\cal{Q}}_1={\rm{ln}}(1 +3\varPhi {\rm{e}}^{-(E_i-\mu_i)/T}+3\bar{\varPhi} {\rm{e}}^{-2(E_i-\mu_i)/T}+{\rm{e}}^{-3(E_i-\mu_i)/T})$, ${\cal{Q}}_2={\rm{ln}}(1+3\bar{\varPhi} {\rm{e}}^{-(E_i+\mu_i)/T} +3\varPhi {\rm{e}}^{-2(E_i+\mu_i)/T}+{\rm{e}}^{-3(E_i+\mu_i)/T})$, and $E_i= \sqrt{\boldsymbol{p}^{\,2}+M_i^2}$ is the dispersion relation. The pressure and energy density can be derived using the thermodynamic relations in the grand canonical ensemble as

For the given temperature $ T $, baryon density $ \rho_{\rm B} $ and asymmetric parameter $ \alpha $, the values of $\phi_{{\rm{u}}},\; \phi_{{\rm{d}}},\; \phi_{{\rm{s}}},\; \varPhi, \;\bar{\varPhi}$, $\mu_{{\rm{u}}}, $ $ \mu_{\rm {\rm{d}}} $ and $ \mu_{{\rm{s}}} $ can be determined by solving the equations derived by minimizing the thermodynamic potential

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

We first present in Fig. 1 the isentropic curves for $\sigma=50,\, 10,\, 5,\,1$ in the $T{\text{-}}\mu_{\rm{B}^{}}$ plane. Each curve in this figure corresponds to an evolutionary trajectory of QGP generated with a specific collision energy. For the two cases of $ \sigma=50 $ and $ 10 $, the fireball only passes through the stable phase during the evolution. However for the cases of $ \sigma=5 $ and 1, both the metastable and unstable phases are involved. The behaviors of σ in the $T {\text{-}} \mu_{\rm{B}}^{}$ plane is mainly determined by the phase structure of the system. In theory, for a fermion system with the temperature and density dependent dynamic mass, there usually exists a first-order phase transition at low temperature. Numerical calculations in different systems show that the qualitative behaviors of $s/\rho_{\rm B}^{}$ are similar to the curves shown in Fig. 1. For the liquid-gas phase transition in nuclear matter, we also derived the similar curves of $s/\rho_{\rm B}^{}$^[48]. Generally, the appearance of the inflection point does not directly corresponds to the position of the phase transition, but the crossover phase transition of QCD at high temperature can induce some perturbation on $s/\rho_{\rm B}^{}$, as shown in the curve of $ \sigma=50 $.

Figure 1.  Isentropic curves for $\sigma=50,\, 10,\, 5,\,1$ for both the symmetric and asymmetric quark matter in the $T {\text{-}} \mu_{\rm{B}}^{}$ plane. (color online)

We plot in Fig. 2 the EOSs of quark matter along different evolutionary paths given in Fig. 1. The numerical results show that for the cases of $ \sigma=1 $ and 5 in symmetric matter and $ \sigma=1 $ in asymmetric matter, the pressure of the system has a maximum and a minimum as the energy density increases, respectively. The region between the maximum and minimum approximately correspond to (lies in exactly) the unstable phase of the system. After passing through the minimum, the pressure grows rapidly with the increase of energy density. For the cases of $ \sigma=5 $ in asymmetric quark matter and $ \sigma=10 $ and 50 in both the symmetric and asymmetric quark matter, the pressure increases monotonically with the rising energy density. Figure 2 also indicates that the equations of state of asymmetric matter are stiffer than those of symmetric matter, although it is not obvious for the case of $ \sigma=50 $. These behaviors of the equations of state under different physical conditions determine the speed of sound in quark matter.

Figure 2.  EOSs along the different evolutionary paths shown in Fig. 1. (color online)

The evolving behaviors of the speed of sound corresponding to the EOSs under different conditions in Fig. 2 are plotted in Fig. 3. This figure demonstrates that the influence of isospin effect on the speed of sound of quark matter. In the chiral restored phase (the high-temperature and/or high-density region), the isospin splitting effect is relatively smaller, because the differences of both the quark masses and chemical potentials between u and d quarks are very small, leading to that the asymmetric parameter $ \alpha $ has a weak effect on the EOSs. The isospin effect is also weak at lower density (lower energy density) in the chiral broken phase with the similar reasons for the chiral restored phase, in particular when the mass of u quark is the same as that of d quark. However, the isospin effect has a relatively larger influence on the EOSs at intermediate density and in the phase transition region. Therefore, from the low-temperature and low-density region to the high-temperature and/or high-density region, the influence of the isospin effect on EOSs goes from weak to strong and then to weak again, which is consistent with the results shown in Fig. 2. With the combination of the above qualitative analysis and the EOSs in Fig. 2, as well as the isetropic trajectories in Figs. 1, 4 and 5, it can be understood that when σ is relatively smaller (e.g., $ \sigma= $1, 5, 10), the speed of sound of asymmetric quark is larger than that of symmetric quark matter in the phase transition region, and a reversal happens when approaching the chiral restored phase. while for $ \sigma=50 $, the sound speed curves of the symmetric and asymmetric quark mater almost coincide because the isentropic curves of them mainly pass through the low-density region and the chiral restored phase at high temperature.

Figure 3.  Square of speed of sound for the EOSs under different conditions along the evolutionary paths shown in Fig. 1. (color online)

To demonstrate the dependence of the speed of sound on temperature and density, we present the contour plots of the speed of sound in both the symmetric quark matter and asymmetric quark matter in Figs. 4 and 5, respectively. From the two figures, it can be seen that in the low-temperature region of the first-order phase transition, the restoration of chiral symmetry causes the square of the speed of sound to quickly approach 1/3, as indicated by the red curve in Fig. 3. In the smooth transition region at high temperatures, the growth of the speed of sound with increasing temperature or density is relatively slower. However, at extremely high temperatures and densities, the square of the speed of sound still approaches 1/3.

The red solid lines in Figs. 4 and 5 correspond to the boundaries where the speed of sound becomes zero with the condition of $ \left(\frac{\partial p}{\partial \epsilon}\right)_{\sigma}=0 $. As pointed out in our previous study^[47], the acoustic wave equation becomes a decaying function within the region enclosed by the red solid lines, and the disturbances cannot propagate outwards within this region. Comparing Figs. 4 and 5, we observe that the region with the negative speed of sound squared in asymmetric matter is smaller than that in symmetric matter, which can also be seen from Fig. 3.

Figure 4.  Contour map of $ c^2_\sigma $ in the $T {\text{-}} \rho_{\rm B}^{}$ plane for symmetric quark matter with $ \alpha=1 $. The red curve is the boundary of vanishing sound speed. (color online)

Figure 5.  Contour map of $ c^2_\sigma $ in the $T {\text{-}} \rho_{\rm B}^{}$ plane for asymmetric quark matter with $ \alpha=1/2 $. The red curve is the boundary of vanishing sound speed. (color online)

Furthermore, combined with the contour plots of the speed of sound under the condition $ \mu_{\rm u}=\mu_{\rm d}=\mu_{\rm{s}} $^[47], it can be observed that the condition $ \rho_{\rm{s}}=0 $ taken in this study has a significant influence on the speed of sound in quark matter, especially in the low temperature and high density region where the chiral symmetry of u and d quarks has restored. The calculation^[47] demonstrates that under the condition of $ \mu_{\rm u}=\mu_{\rm d}=\mu_{\rm{s}} $, the speed of sound at low temperatures rapidly increases to nearly $ \sqrt{1/3} $ after the restoration of chiral symmetry of u and d quarks. However, as the density increases, the presence of strange quarks causes the equation of state to become softer, resulting in a subsequent decrease of the sound speed, which tends to be $ \sqrt{1/3} $ again after the restoration of chiral symmetry of s quark at higher densities.

To compare the value of the speed of sound in symmetric and asymmetric quark matter, we plotted the contour maps of $ r_{\rm{s}}=c^2_{\sigma=1/2} $/$ c^2_{\sigma=1} $ in Fig. 6. In the gray area, the speed of sound in asymmetric quark matter is slightly larger than that in symmetric quark matter, and $ r_{\rm{s}}=1 $ on the boundary of the gray area. In the red region, $ r_{\rm{s}}<1 $, it means that the value of the speed of sound in asymmetric quark matter is smaller than that in symmetric quark matter. It should be noted that the difference in speed of sound between the two cases in the gray area is very small ($ r_{\rm{s}}<1.005 $ in most of the gray region). The negative ratio appears in the region between the two boundaries (black solid curve and dashed curve) of the vanishing sound speed. In this region, the square of the speed of sound in asymmetric quark matter is positive, while in symmetric quark matter, the mechanical instability violates the sound wave equation with the $ c^2_{\sigma} $ being negative. The region where the ratio takes a larger absolute value locates near the boundary of (black solid curve) of the vanishing sound speed in symmetric quark matter.

Figure 6.  Contour map of $r_{\rm{s}}=c^2_{\sigma=1/2}$/ $ c^2_{\sigma=1} $ in the $T {\text{-}} \rho_{\rm B}^{}$ plane. The solid (dashed) black curve is the boundary of vanishing sound speed in symmetric (asymmetric) matter. In the gray area, the squared speed of sound in asymmetric matter is slightly larger than that in the symmetric matter with $r_{\rm{s}} > 1$; $r_{\rm{s}}=1$ on the boundary of the gray area; $r_{\rm{s}} < 1$ in the red area. (color online)

We derived the formula for calculating the speed of sound in asymmetric quark matter in the temperature and density space, and calculated the speed of sound in the PNJL quark model. The results indicate that the boundary of vanishing sound speed in asymmetric quark matter is smaller than that in symmetric quark matter, meaning that the range where the acoustic wave equation fails in asymmetric quark matter is smaller than that in symmetric quark matter. The calculations also indicate that in most of the stable phase, the speed of sound in asymmetric quark matter is slightly larger than that in symmetric quark matter. Additionally, it should be pointed out that only the change of sound speed caused by isospin asymmetry of u, d quark chemical potential is considered in this study. The difference in u, d quark dynamic masses is not included. After introducing the isospin interactions that cause the splitting of u, d quark masses, it is expected that the isospin effects on the speed of sound will be more pronounced. The relevant research is currently underway.

Acknowledgments The work was supported by National Natural Science Foundation of China (Grant No. 11875213) and Natural Science Basic Research Plan in Shaanxi Province of China (Program No. 2024JC-YBMS-018).