Charmed meson structure across crossover from SU(4) Polyakov quark meson model with isospin asymmetry

Recent progress in ultra-relativistic heavy-ion experiments–such as the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory (BNL), Large Hadron Collider (LHC) at CERN, upcoming Compressed Baryonic Matter (CBM) experiment at the Facility for Antiproton and Ion Research (FAIR), and Nuclotron-based Ion Collider fAcility (NICA) at the Joint Institute for Nuclear Research (JINR) –has significantly advanced the investigation of strongly interacting matter. Quantum Chromodynamics (QCD), the fundamental theory of strong interactions among quarks and gluons, predicts that at high temperatures and baryon densities, hadronic matter undergoes a phase transition into a deconfined and chirally symmetric phase known as the quark-gluon plasma (QGP) [1–6].

In recent decades, effective low-energy QCD models have been developed to investigate different features of the QCD phase transition. These approaches incorporate the key concepts of spontaneous chiral symmetry breaking and color-charge deconfinement [7–9]. Chiral symmetry is a fundamental symmetry in QCD in the limit of vanishing quark masses [7, 8]. For finite quark masses, chiral symmetry becomes spontaneously broken, with the chiral condensate serving as the associated order parameter. When quark masses are small (up u- and down d-quarks), chiral symmetry becomes explicitly broken, resulting in the pseudo-Goldstone bosons acquiring finite masses; for example, the light pseudoscalar mesons as pions and kaons. For heavier quarks, such as strange and charm quarks, the explicit breaking is stronger, leading to the formation of hidden and open charmed mesons, such as D-resonances and $ \chi_{c0} $. At high temperature or density, the chiral condensate melts, signaling the restoration of chiral symmetry and the onset of meson state degeneracy.

Low-energy QCD approaches have been widely employed to study the QCD phase structure. These approaches, such as chiral perturbation theory [10], Dyson-Schwinger Equations (DSEs) [11–13], hadron resonance gas, and other statistical models, such as [14–16], the Polyakov Nambu-Jona-Lasinio (PNJL) model [17–20], Polyakov linear sigma model (PLSM), or Polyakov quark meson (PQM) model [21–24], have been widely applied. In the world of $ N_f $ quark flavors with u- and d-quarks, i.e., SU(2) flavor symmetry [25–27], pions are formed, and the light quark condensate is investigated. The extension to SU(3) symmetry includes the strange s-quark, allowing for the investigation of the strange condensate, nonet meson masses, and QCD phase transition [28–30]. Further extension to SU(4) flavor symmetry incorporates the charm quark and its associated condensate, making it possible to analyze the charmed meson states in a thermal and dense QCD medium [31–34]. The validity of the PQM model has been confirmed by extensive comparisons with lattice QCD simulations, demonstrating key features of QCD thermodynamics and phase transitions [35–37], including the non-extensive effects on the QCD equation of state (EoS) [38], SU(3) isospin asymmetry [39], SU(4) symmetry and QCD thermodynamics [33, 34], and viscous properties [36]. Further extensions studied the chiral symmetry containing the vector meson states at finite density [22] and magnetic field [40]. Other studies explored the QCD EoS in the presence of a magnetic field [41–43] and magnetic catalysis [44, 45].

One of the main targets of the research programs at NICA (JINR) [46] is to investigate in-medium properties of the QCD phase diagram at high density and moderate temperatures by colliding heavy ions, exploring the onset of deconfinement, locating the chiral symmetry restoration boundary, and searching for the critical endpoint (CEP). The created matter will naturally exhibit finite isospin asymmetry, introducing a nonzero isospin chemical potential ($ \mu_I $) into the system. Theoretical approaches, such as the PQM model and lattice QCD, provide essential guidance for analyzing experimental data, accurately mapping the chiral QCD phase structure, and identifying the CEP in iso-asymmetric QCD matter. The motivation behind the current manuscript is to utilize and examine the inclusion of finite isospin asymmetry and the charm chiral condensate on the chiral QCD phase structure, the chiral phase diagram of QCD matter, thermodynamic observables, diagonal and off-diagonal fluctuations of baryon-charm interactions, and the mass spectrum of (pseudo)scalar meson states at finite temperature T, particularly at very high temperatures on the order of $ 1\; \text{GeV} $. The inclusion of the charm chiral condensate significantly improves the capability of the PQM model to evaluate the QCD phase structure, especially at high temperatures. This enhancement is relevant not only for particle physics but also for modeling the EoS in hyperonic matter, Bose-Einstein condensation (BEC), color superconductivity (CSC), and EoS for the properties of neutron stars (NS) and merger remnants [47, 48].

The remainder of this manuscript is organized as follows. Section II outlines the structure of the SU(4) PQM model in finite isospin asymmetry. Section III focuses on the analysis of the SU(4) PQM model in a thermal and dense QCD medium. The PQM order parameters and chiral QCD phase diagram are discussed in subsection III.A, while thermodynamic observables, the EoS of QCD matter, and fluctuations and correlations of conserved charges over a wide temperature range are explored in subsection III.B. Subsection III.C examines the meson mass spectrum at both zero and finite temperatures. Finally, Section IV provides a summary and conclusion.

The SU($ N_f $) PQM model addresses the phase transition between the confinement and deconfinement phases. The chiral Lagrangian is incorporated with the Polyakov-loop potential as $ {\cal{L}}_{{\rm{PQM}}} = {\cal{L}}_{{\rm{chiral}}} - {\cal{U}}(\phi, \bar{\phi}, T) $. The chiral Lagrangian $ {\cal{L}}_{{\rm{chiral}}} $ for ($ N_f=4 $) quark flavors and ($ N_c = 3 $) color degree of freedom is given as follows:

where

Here, the first term on the right-hand side (rhs) of Eq. (1) represents the Lagrangian density for the baryonic (fermionic) field $ \psi_f $, where $ D_{\nu} $, ν, $ \gamma^{\nu} $, and g denote the covariant derivative, Lorentz index, chiral spinors, and Yukawa coupling constant, respectively. The index f refers to the quark flavors, given as $ f = [(l = u, d), s, c] $. The remaining terms on the rhs of Eq. (1) correspond to the contributions from the mesonic (bosonic) field $ {\cal{M}} $. The generator operator $ {\cal{M}} $ is represented by a complex matrix [49, 50] as

Within the U(4) algebra framework with the pseudoscalar ($ \pi_a, J^{{\cal{PC}}=0^{-+}} $) and scalar ($ \sigma_a, J^{{\cal{PC}}=0^{++}} $) sectors, the generators $ T_a= \hat{\lambda}_a/2 $ can be derived from Gell-Mann matrices $ \hat{\lambda}_a $ [52], where the index $ a=0,\cdots,\,15 $ (see App. A). The fields $ \sigma_a,\, \pi_a $ denote the scalar ($ J^{PC} = 0^{++} $) and pseudoscalar ($ J^{PC} = 0^{-+} $) nonet meson fields, respectively [49, 50]. For the sake of clarity, Eq. (3) illustrates the meson sectors given in Table 1. In vacuum phenomenology, the algebraic group $ U(4) $ symmetry is decomposed as $ U(4) \approx SU(4) \times U(1)_A $, where $ U(1)_A $ is represented by the anomaly term $ {\cal{L}}_{U(1)_A} = {\cal{C}}({\rm{Det}}[{\cal{M}}] + {\rm{Det}}[{\cal{M}}^{\dagger}]) $. Moreover, one can assume that $ H \neq 0 $, $ {\cal{C}} \neq 0 $, and $ \lambda \neq 0 $. This implies that the chiral Lagrangian leads to nonzero quark masses, and the anomaly term is explicitly broken. As a consequence of chiral symmetry spontaneously breaking, the mean value of the mesonic field, $ \langle {\cal{M}} \rangle = T_0 \sigma_0 + T_3 \sigma_3 + T_8 \sigma_8 + T_{15} \sigma_{15} $, must carry specific quantum numbers, with only $ h_0 $, $ h_3 $, $ h_8 $, and $ h_{15} $ being non-vanishing. To address finite isospin asymmetry, where the quark flavor masses are not completely degenerate (i.e., $ m_u \neq m_d,\;m_s,\;m_c $), we utilize an orthogonal basis transformation. This transformation converts the condensates from the original basis, $ \sigma_0 $, $ \sigma_3 $, $ \sigma_8 $, and $ \sigma_{15} $, to the pure up ($ \sigma_u $), down ($ \sigma_d $), strange ($ \sigma_s $), and charm ($ \sigma_c $) flavor basis, respectively. This conversion can be expressed as

Pseudoscalar	–pion triplet states $ \pi^\pm=(\pi_1\pm {\rm i} \pi_2)/\sqrt{2},\pi^0=\pi_3 $,
sectors	–four kaon states $ K^\pm=(\pi_4\pm {\rm i} \pi_5)/\sqrt{2}, K^0 =(\pi_6+{\rm i} \pi_7)/\sqrt{2}\quad $ and $ \overline{K}^0=(\pi_6 -{\rm i} \pi_7)/\sqrt{2} $,
	–isoscalar fields $ \eta_{\ell,s} $ mix and generate the physical fields η and $ \eta^\prime $ [21, 30, 51].
$ J^{PC} = 0^{-+} $	–open charmed D resonance $ D^\pm =(\pi_{11} \pm i \pi_{12})/\sqrt{2}, D^0 =(\pi_{9} - {\rm i} \pi_{10})/\sqrt{2}, \overline{D}^0=(\pi_{9} +{\rm i} \pi_{10})/\sqrt{2} $,
	–open strange-charmed $ D_S^\pm=(\pi_{13} \pm {\rm i} \pi_{14})/\sqrt{2} $, and
	–hidden charmed $ \eta_C=(\pi_0-\sqrt{3} \pi_{15})/2 $.
Scalar	–isotriplet of $ \vec{a}_0 $ meson $ a_0^\pm=(\sigma_1 \pm i \sigma_2)/\sqrt{2}, a_0^0=\sigma_3 $
sectors	–four kaon states $ \kappa ^\pm=(\sigma_4 \pm{\rm i} \sigma_5)/\sqrt{2} $, $ \kappa^0 =(\sigma_6+{\rm i}\sigma_7)/\sqrt{2} $ and $ \overline{\kappa}^0=(\sigma_6 - {\rm i} \sigma_7)/\sqrt{2} $,
	–isoscalar fields $ \sigma_{\ell,s} $ mix and generate the physical fields $ f_0 $ and σ [21, 30, 51],
$ J^{PC} = 0^{++} $	–open charmed $ D_0 $ resonance $ D_0^\pm =(\sigma_{11} \pm {\rm i} \sigma_{12})/\sqrt{2}, D_0^0 =(\sigma_{9} -{\rm i} \sigma_{10})/\sqrt{2},\overline{D}_0^0=(\sigma_{9} + {\rm i} \sigma_{10})/\sqrt{2} $,
	–open strange-charmed $ D_{S,0}^\pm=(\sigma_{13} \pm {\rm i} \sigma_{14})/\sqrt{2} $, and
	–hidden charmed $ \chi_{c0}=(\sigma_3-\sqrt{3} \sigma_{15})/2 $.

Table 1.  Summary of the field contents of meson sectors.

The chiral symmetry is explicitly broken by a $ 4\times 4 $ matrix, $ H= T_a h_a $, where the explicit symmetry breaking terms, $ h_a $,

can be determined from the partially conserved axial current (PCAC) relations [39, 49]. The parameters of explicit symmetry breaking for light ($ l=u\text{-}, d $-quark), strange (s-quark) and charm (c-quark) are given as

In this regard, we utilize the PQM model to analyze the QCD phase structure in finite isospin asymmetry. First, we investigate the impacts of isospin asymmetry to distinguish between the light condensates of u- and d-quarks. This arises due to the breaking of isospin asymmetry, where $ \sigma_3 = a_0^0= \big(f_{K^\pm}-f_{K^0}\big)\neq 0 $ [39]. Second, we define the explicit symmetry breaking for light quarks $ (l=u,d) $ as $ h_u = h_l+h_3 $ and $ h_d = h_l-h_3 $, where the term $ h_3 $ representing explicit symmetry breaking is given by $ h_3 = m_{a0}^2 \big(f_{K^\pm}-f_{K^0}\big) $ [39]. Finally, in the term $ -2{\rm{Tr}} [\varepsilon ({\cal{M}}^{\dagger} {\cal{M}})] $ in Eq. (1), ε represents the direct contribution of bare quark to nonet meson masses:

For the sake of simplicity, one can set $ \varepsilon_u = \varepsilon_d = \varepsilon_s = 0 $ [30, 31]. The used parameters $ m^2 $, $ h_u $, $ h_d $, $ h_s $, $ h_c $, $ \lambda_1 $, $ \lambda_2 $, and $ {\cal{C}} $ are expressed in terms of $ m_\sigma $ [21]. In the current manuscript, the relevant parameters assume $ m_\sigma=800 $ MeV. [21]. The Polyakov-loop potential, $ {\cal{U}} (\phi, \overline{\phi}, T) $, incorporates degrees of freedom of color-charges and the dynamics of the color gluon-quark interaction into the chiral Lagrangian of the PQM model. The potential is estimated by calculating the thermal expectation value of a color–traced Wilson loop in temporal space [53, 54] as $ \phi=\rm{Tr}_c[{\cal{L}}]/N_c $ and its Hermitian conjugate $ \overline{\phi}=\rm{Tr}_c[L^\dagger]/N_c $, where $ {\cal{L}} $ and $ {\cal{L}}^\dagger $ are matrices in the representation of the SU($ N_c $) color-charge gauge group [24, 55, 56]:

where $ A_0 $ and $ {\cal{P}} $ are the temporal component and path ordering of the gluon gauge field, respectively. Various functional forms of the Polyakov loop potential exist [21, 57–60]; in this work, we adopt the Polyakov loop potential proposed by Fukushima [61], inspired by a strong-coupling interaction

where $ a=664.0\; $MeV and $ b=196.2\; $MeV are fixed parameters. The parameter a controls the confinement transition in pure gauge, while b is related to the mixing of the chiral and deconfinement transition [61]. The expectation value of the Polyakov loops is finite and approaches unity (indicating the deconfined phase) at high temperatures, while it diminishes to zero or takes small values at low temperatures (indicating the confined phase) [62]. Hence, this order parameter serves as a reliable indicator of the confinement-deconfinement phase transition.

The grand canonical partition function $ {\cal{Z}} $, which depends on the temperature T and chemical potentials of f-th quark flavor $ \mu_f $, is defined by the path integral over all fermions, anti-fermions, and bosons. In mean-field approximation [63], all fields are treated as averages in space ($ \vec{x} $) and imaginary time (τ) as

where $ {\cal{N}}_f ~ \rm{and} ~ \hat{{\cal{H}}} $ are defined by the quark number density and the chiral Hamiltonian density, respectively. The chemical potentials $ \mu_f $ are related to the conserved quantum baryon (B), electric charge (Q), isospin (I), strangeness (S) and charm (C) numbers for each of the quark flavors:

The exchanges between particles and antiparticles are incorporated through the grand canonical partition function ($ {\cal{Z}} $), which is constructed from the free energy density as

The first term on the rhs of Eq. (11) represents the tree-level mesonic potential. It accounts for the contribution of the pure multiplet of the nonet meson states with a basis of light ($ \sigma_u $, $ \sigma_d $), strange ($ \sigma_s $), and charm ($ \sigma_c $) condensates. The purely bosonic (mesonic) potential is given as

The second term on the rhs of Eq. (11) represents the fermion contributions to the PQM model [22, 61] as

where the first term on the rhs of Eq. (13) corresponds to the fermionic vacuum contribution; this term is usually neglected in calculations. Moreover, one can regularize the dimensions to avoid the diverging integral [64, 65]:

where Λ is defined as a regularisation scale parameter. Meanwhile, the second term on the rhs of Eq. (13) gives the fermionic contribution in a thermal and dense QCD medium, whereas $ g_D^{(+)}(T,; \mu_f) $ and $ g_D^{(-)}(T,; \mu_f) $ are defined as the Fermi-Dirac distribution functions in the presence of expectation values of Polyakov loop variables, following widely used notations:

where $ E_f^{(\pm)}=E_f(k) \mp \mu_f $ represent the energy-momentum dispersion relations, with the upper sign applying to quarks and the lower sign to antiquarks. The energy-momentum dispersion relations for quark flavors are defined as $ E_f(k)=(k^2+m_f^2)^{1/2} $, where $ m_f $ denotes the constituent quark masses:

The term $ g_D^{(-)}(T,; \mu_f) $ can be evaluated by substituting $ E_f^{(+)} $ with $ E_f^{(-)} $ and replacing the Polyakov loop variable ϕ with its conjugate $ \overline{\phi} $, and vice versa. As the chiral symmetry is spontaneously broken, the PQM Lagrangian undergoes a shift in the meson fields $ \sigma_f $ by their vacuum expectation values, denoted as $ \sigma_f \longrightarrow \sigma_f + \langle \sigma_{f0} \rangle $. The values of quark chiral condensates in a vacuum are calculated in terms of decay widths of pion, kaon, and D-meson using the partially conserved axial-vector current (PCAC) relations $ \overline{\sigma}_{\ell o} = f_\pi = 92.4\; $MeV, $ \overline{\sigma}_{so} =(2 f_K - f_\pi)/ \sqrt{2}= 94.5\; $MeV, and $ \overline{\sigma}_{co} =(2 f_D - f_\pi)/\sqrt{2} = 293.87\; $MeV. We use the experimental values [66, 67], where $ f_{\pi^{\pm}} = f_{\pi^0} = 92.4\; {\rm{MeV}} $, $ f_{K^{\pm}} = 113\; {\rm{MeV}}, ~ f_{K^0} =113.453 \;{\rm{MeV}} $ and $ f_D= 254\; $MeV are the decay constants of pion, kaon, and D meson states, respectively. Within the mean-field approximation (MFA), one can examine the thermodynamic observables that characterize the chiral QCD phase structure of quark-hadron matter as functions of independent parameters, such as temperature, chemical potentials, and magnetic field. These observables are derived by examining the thermal behavior of the order parameters of the PQM model. Within this framework, the quark chiral condensates $ \overline{\sigma}_u $, $ \overline{\sigma}_d $, $ \overline{\sigma}_s $, and $ \overline{\sigma}_c $, as well as the expectation values of the Polyakov-loop variables $ \overline{\phi} $ and $ \overline{\bar{\phi}} $, are determined by locating the global minimum of the real part of the PQM thermodynamic free energy density $ {\cal{R}} \;[{\cal{F}}_{\rm PQM}(T, \mu_f)] $, as given by Eq. (11), to the relevant order parameters at the saddle point,

The Polyakov loop is a complex-valued matrix that renders the PQM free energy density in Eq. (11) a functional of complex variables. Notably, the physical observables derived from the imaginary part of this function vanish, as discussed in Ref. [68]. This involves solving a set of equations simultaneously. The thermal evolution of the order parameters is analyzed using the gap equations and by tracking the approach towards the saddle point.

This section investigates the properties of the quark-hadron phase structure using the SU(4) PQM model. The analysis highlights several thermodynamic observables, such as chiral condensates, deconfinement order parameters, pseudo-critical temperatures, and other thermodynamic quantities, compared with lattice QCD simulations. Additionally, we examine the quark susceptibilities, fluctuations, and correlations of conserved charges associated with the charm quantum number. Finally, the mass spectrum of nonet meson states is introduced within the framework of the SU(4) PQM model.

It is helpful to present several calculations to illustrate the impact of isospin asymmetry on the chiral phase transition and deconfinement order parameters in the PQM model for $ N_f = 2+1+1 $ quark flavors. To achieve a reliable differentiation between u- and d-quark condensates, it is necessary to estimate the influence of finite isospin asymmetry on the SU(4) PQM order parameters. To this end, we begin with the real part of the PQM thermodynamic free energy density in the mean-field approximation (MFA) (Eq. (11)). This term should be globally minimized at a saddle point with respect to the corresponding mean field, as outlined in the gap equations (Eq. (17)). The solutions of the gap equations provide a complete set of equations that illustrate the evolution of the chiral condensates $ \sigma_f $ with $ [f=(l=u,d),s,c] $ quark flavors, as well as the order parameters of Polyakov loop potential ϕ and $ \bar{\phi} $, in the thermal and dense QCD medium.

Figure 1 depicts the quark chiral condensates $ \sigma_f/\sigma_{f_0} $ normalized to their values in a vacuum and the order parameters of the Polyakov loop variables (ϕ and $ \bar{\phi} $) as a function of temperature at different values of isospin ($ \mu_I $) and baryon ($ \mu_B $) chemical potentials: left panel (a) at $ \mu_I = 0.0 $, middle panel (b) at $ \mu_I = 150.0\; \text{MeV} $, and right panel (c) at $ \mu_I = 225.0\; \text{MeV} $. In the current calculations, we assume that $ \mu_B = 2.0\mu_I $, and other chemical potentials for conserved charges are set to zero. The left panel (a) displays the normalized expectation values of the nonstrange ($ \sigma_{\ell} $, solid curve), strange ($ \sigma_s $, dashed curve), and charm ($ \sigma_c $, dotted curve) chiral condensates, as well as the order parameter for the Polyakov loop variables $ \phi = \bar{\phi} $ (dotted-dashed curve), at zero isospin chemical potential ($ \mu_I = 0.0 $). The nonstrange condensate is defined as $ \sigma_\ell = (\sigma_u + \sigma_d)/2 $. In the hadronic phase, when the temperature is below the pseudo-critical temperature $ T_\chi $, the values of $ \sigma_{f}/\sigma_{f_0} $ start at unity and remain unchanged until the phase transition region around $ T \approx T_\chi $. This indicates that the expectation values of the quark chiral condensates are nearly indistinguishable, signifying a confined state. As the system's temperature increases, the chiral condensates gradually decrease, resulting in distinct curves that characterize the nonstrange, strange, and charm chiral condensates. The smooth decrease indicates a slow transition, or ''crossover''. At high temperatures, the chiral condensates are significantly suppressed, indicating a transition to a deconfined phase in the QCD system. The thermal evolution of the deconfinement phase transition is evident in the behavior of the Polyakov loop order parameters ($ \phi,\; \bar{\phi} $). The order parameters of the Polyakov loop variables are identical, i.e., $ \phi = \bar{\phi} $, at $ \mu_f=0 $. At low temperatures (confined phase), the confined hadron phase dominates, and the order parameters of the Polyakov loop variables vanish or take small values. As the temperature increases, the confined phase transforms smoothly to the deconfined phase, occurring at a deconfinement temperature $ T_\phi $. At higher temperatures (deconfined phase), the order parameters of the Polyakov loop variables for the Fukushima-type potential approach unity, which is related to the gluonic dynamics and static color-charge free energy.

Figure 1.  (color online) Left panel (a) illustrate the thermal behaviors of the PQM order parameters ($ \sigma_f/\sigma_{f_0} $, ϕ, and $ \bar{\phi} $) at $ \mu_I=0.0 $. The middle and right panels (b) and (c) show the same as (a) but at $ \mu_I=150.0\; $MeV and $ \mu_I=225.0\; $MeV, respectively.

In the middle (b) and right (c) panels of Fig. 1 at finite $ \mu_I $, the significance of the isospin parameters ($ \sigma_3 $, $ h_3 $) becomes evident in distinguishing between the $ u $- (dash-double-dotted curves) and $ d $- (long-dashed-dotted curves) quark condensates. During the phase transition, the normalized light quark condensates split into two separate curves representing the u- and d-quark condensates. It is apparent that the pseudo-critical temperature $ T_\chi $ decreases as $ \mu_I $ increases. Additionally, the normalized strange and charm condensates remain unaffected by changes in $ \mu_I $, and the value of $ T_\chi $ decreases as $ \mu_B $ increases. Furthermore, the nonstrange and strange chiral condensates exhibit greater sensitivity to $ \mu_I $ and $ \mu_B $ than the charm condensate. It is noteworthy that both the isospin $ \mu_I $ and baryon $ \mu_B $ chemical potentials lead to the same conclusion: the phase transition occurs at a lower temperature $ T_\chi $. When the isospin chemical potential $ \mu_I $ is introduced, the thermal evolution of the order parameters associated with the deconfinement transition becomes more pronounced. Specifically, the Polyakov loop variable ϕ (long-dashed curves) and its conjugate $ \bar{\phi} $ (triple-dashed curves) begin to differ, i.e., $ \phi \neq \bar{\phi} $. It is evident that the thermal behavior of ϕ consistently exceeds that of $ \bar{\phi} $ as $ \mu_I $ increases. Increasing either $ \mu_I $ or $ \mu_B $ results in a gradual rise in the Polyakov loop variables and thereby enhances the deconfinement transition. It is significant that the effects of both $ \mu_I $ and $ \mu_B $ lead to the same outcome: the deconfinement transition occurs at a lower temperature $ T_\phi $.

Figure 2 illustrates the behavior of the normalized chiral condensates ($ \sigma_f/\sigma_{f_0} $) and the order parameters of the Polyakov loop variables (ϕ and $ \bar{\phi} $) as functions of the isospin chemical potential ($ \mu_I $) at different temperatures: $ T = 120.0\; \text{MeV} $ (left panel), $ T = 145.0\; \text{MeV} $ (middle panel), and $ T = 165.0\; \text{MeV} $ (right panel). The impact of temperature (T) and isospin chemical potential ($ \mu_I $) on the light quark condensates ($ \sigma_u $ and $ \sigma_d $) and heavy quark condensates ($ \sigma_s $ and $ \sigma_c $) highlights significant differences between the light and heavy quark sectors. The chiral condensates ($ \sigma_u $, $ \sigma_d $, and $ \sigma_s $) show pronounced sensitivity to both T and $ \mu_I $. The introduction of isospin asymmetry causes the chiral condensates of the light quark sector to become distinguishable: as $ \mu_I $ increases, $ \sigma_u $ decreases more rapidly than $ \sigma_d $ and $ \sigma_s $. In the beginning, the chiral condensates remain at their vacuum values ($ \sigma_{f_0} $) at low $ \mu_I $, where the system resides in a confined phase with chiral symmetry spontaneously broken. As $ \mu_I $ increases, the chiral condensates of the $ \ell $- and s-quark sectors gradually decrease due to the restoration of the chiral symmetry, with a suppression near the phase transition. However, the thermal behavior of the chiral condensates exhibits a smooth phase transition or smooth crossover. This suppression becomes more pronounced with increasing T. While the charm quark condensate ($ \sigma_c $) is not sensitive to the isospin chemical potential, $ \sigma_c $ remains nearly fixed as both T and $ \mu_I $ increase.

Figure 2.  (color online) The same as Fig. 1 but varying the isospin chemical potential at different values of temperature T = 120.0 MeV (left panel), 145.0 MeV (middle panel), and 165.0 MeV (right panel).

At lower temperatures, below the critical endpoint (CEP), i.e., $ T = 120.0\; \text{MeV} $ (left panel), the system remains in the confined phase, characterized by high-value chiral condensates and small values of the order parameters of the Polyakov loop variables. As $ \mu_I $ increases, the transition towards the deconfined phase becomes more pronounced, with a sharper decline in chiral condensates and a significant rise in the Polyakov loop order parameters. At temperatures above the CEP, i.e., $ T = 165.0\; \text{MeV} $ (right panel), the system is almost entirely in the deconfined phase. Under $ \mu_I $ dependence, the chiral condensates are suppressed, while the Polyakov loop order parameters are elevated. The behavior of the order parameters of the Polyakov loop variables (ϕ and $ \bar{\phi} $) also exhibits isospin asymmetry. At low $ \mu_I $, the two parameters are nearly equal, reflecting a symmetric QCD medium. However, as $ \mu_I $ increases, ϕ rises more rapidly than $ \bar{\phi} $. Furthermore, increasing T shifts the deconfinement transition to higher values, indicating that the system transitions more rapidly from the confined to deconfined phase. The thermal behavior of ϕ is more pronounced than that of its conjugate, $ \bar{\phi} $. When $ \mu_I $ is introduced, the evolution of the order parameters for the deconfined phase transition becomes distinguishable. It is important to highlight that both temperature (T) and isospin chemical potential ($ \mu_I $) are similar, resulting in a lowered deconfinement transition temperature, $ T_\phi $. The pseudo-critical temperature ($ T_{\chi} $) at different values of the isospin chemical potential can be estimated by identifying the location of the inflection point as follows:

The obtained results map out the chiral phase diagram ($ T/T_{\chi}, \mu_I/m_\pi $) plane at $ \mu_B = 0\; $MeV. Figure 3 illustrates the dependence of the normalized pseudo-critical temperature ($ T/T_{\chi} $) on the normalized isospin chemical potential to the pion mass ($ \mu_I/m_\pi $) for the SU(4) PQM model. The current result (solid black curve) is compared with lattice QCD calculations from Refs. [69] (shaded band) and [70] (open symbols). The coordinates of ($ T/T_{\chi}, \mu_I/m_\pi $) were estimated by identifying the temperature corresponding to the inflection point of the chiral susceptibility at vanishing chemical potentials. This result reveals that the chiral phase diagram exhibits a decreasing trend in $ T_{\chi} $ with increasing $ \mu_I $. Furthermore, the obtained PQM curve shows good agreement with lattice QCD calculations [69, 70], particularly for $ \mu_I/m_{\pi} < 1 $. At higher values of $ \mu_I $, where lattice QCD data are absent, the PQM model provides an extrapolation, predicting a continued suppression of the pseudo-critical temperature as $ \mu_I $ increases. The PQM curve smoothly extends towards lower temperatures and higher isospin chemical potentials. This agreement validates the effectiveness of the PQM model in capturing essential QCD thermodynamic observables over a wide range of isospin chemical potentials. The SU(4) PQM model and lattice QCD simulations have demonstrated that increasing the isospin chemical potential ($ \mu_I $) leads to a decrease in the pseudo-critical temperature ($ T_\chi $) associated with chiral symmetry restoration. This finding suggests that chiral symmetry is restored at lower temperatures under such conditions, potentially leading to a shift in the location of the critical-endpoint (CEP) within the QCD phase diagram. In heavy-ion collisions at NICA, the produced matter is expected to exhibit finite isospin asymmetry due to the nature of of the colliding heavy ions, resulting in a nonzero isospin chemical potential $ \mu_I $.

Figure 3.  (color online) Chiral phase diagram ($ T-\mu_I/\mu_\pi $) at vanishing baryon chemical potential. The PQM results (solid curves) are confronted with lattice QCD calculation [69] (shaded band) and [70] (open symbols).

Thermodynamic observables play an essential role in characterizing of the chiral–deconfinement phase transition for $ N_f = 2+1+1 $ quark flavors in a thermal and dense QCD medium. For a statistical system in equilibrium within volume V, the thermodynamic pressure can be estimated from the grand canonical free energy density as $ p(T, \mu_f) = {\cal{F}}_{{\rm{PQM}}}(T, \mu_f) $, and the normalized interaction measure (trace anomaly) is given by

where the energy density is given by $ \epsilon (T,\,\mu_f)= T^2 \partial \big[p/T\big]/ \partial T = -p(T,\,\mu_f) + T\,s(T,\,\mu_f) $, while the entropy density is defined as $ s(T,\,\mu_f)=\partial p/\partial T $. In the ideal gas limit, the Stefan–Boltzmann (SB) limits can be determined from the grand canonical partition function for $ (N_c^2 - 1) $ massless gluons and $ N_f $ massless quark flavors, using lowest-order perturbation theory [71]:

where the first term on the rhs of Eq. (20) represents the gluonic contribution, the second term corresponds to the fermionic contribution, and the remaining terms account for the contributions of an ideal gas at finite chemical potential. Figure 4 presents various dimensionless thermodynamic observables: the normalized pressure $ p/T^4 $, interaction measure (trace anomaly) $ I/T^4 $, and entropy density $ s/T^3 $ (left panel (a)), along with the quark flavor susceptibilities $ \chi_f^2 $ (right panel (b)) as functions of $ T/T_\chi $ at $ \mu_f = 0.0\; \text{MeV} $. The results are obtained from the SU(4) PQM model in the presence of the Polyakov loop variables (solid curves). The PQM results show good agreement with lattice QCD simulations [72–74].

Figure 4.  (color online) A comparison between the normalized PQM thermodynamic observables $ p/T^4,\, I/T^4 $, and $ s/T^3 $ (solid curves) at $ \mu_f=0.0$ MeV for SU(4) quark flavors with the lattice QCD simulations [72] (closed circles) and [73] (open circles) (left panel (a)). The quark flavor susceptibilities $ \chi_f^2 $ are confronted with the corresponding lattice QCD calculations [74] (closed symbols) (right panel (b)).

At low temperatures, the normalized pressure is suppressed due to confinement and then gradually increases with $ T/T_\chi $ as new degrees of freedom emerge during the transition to the deconfined phase. At high temperatures, the pressure closes to the SB limit, indicating that the system behaves like an massless ideal gas of free color quarks and gluons, where the interactions are effectively negligible.

The increase in energy density with rising $ T/T_\chi $ reflects the rapid liberation of degrees of freedom during the QCD phase transition, while the entropy density quantifies the disorder of the system during the hadron-quark transition. Furthermore, the interaction measure (trace anomaly) quantifies deviations from conformal symmetry, where the energy-momentum tensor $ T^{\mu\nu} $ vanishes in the limit of a massless ideal gas. The trace anomaly shows a pronounced peak around $ T/T_\chi \approx 1 $, marking the hadron-quark phase transition. This peak corresponds to the maximal interaction effects during the crossover to the quark-gluon plasma (QGP) phase. At high temperatures, the trace anomaly decreases and approaches zero, indicating the recovery of approximate conformal symmetry.

The grand canonical ensemble treats conserved charges as external input parameters. During a heavy-ion collision, the baryon number, electric charge, and strangeness quantum numbers remain conserved. In QCD, conserved charges are associated with four quark chemical potentials that can be introduced in a $ 2+1+1 $ quark flavor system. The set of chemical potentials $ \mu_f $ is defined in Eqs. (10). The higher-order susceptibilities (or fluctuations) with respect to the quark chemical potentials can be expressed in terms of derivatives of the normalized pressure in the PQM model with respect to the corresponding quark chemical potential:

where $ \hat{\mu}_f \equiv \mu_f/T $, ensuring that the cumulants remain dimensionless, and the indices i, j, and k denote the order of partial derivatives. Following the standard procedure in lattice QCD calculations, various fluctuations (diagonal) and correlations (off-diagonal) of different conserved quantum numbers, $ x = [B,\, Q,\, S,\, C] $, can be estimated from the normalized pressure in the PQM model with respect to the corresponding chemical potential of conserved quantum number:

The superscripts α, β, γ, and δ represent the orders of the derivatives, where $ \hat{\mu}_x \equiv \mu_x/T $. Equations (21) and (22) illustrate the (diagonal) fluctuations and (off-diagonal) correlations of quark flavors and conserved quantum numbers, respectively. The moment products or diagonal (off-diagonal) susceptibilities can be established by taking derivatives of the normalized thermodynamic pressure of the SU(4) PQM model with respect to the corresponding dimensionless chemical potentials, as defined in Eq. (10). These moment products can be related to experimentally measured particle multiplicities, encompassing the mean, variance, skewness, kurtosis, and higher-order cumulants, as discussed in [75].

In this section, we present the quark flavor susceptibilities $ \chi_f^2 $ with respect to the light $ \ell = \{u,\,d\} $, strange s, and charm c quark flavors as functions of temperature T at vanishing $ \mu_f $. Furthermore, we show the thermal behavior of diagonal (off-diagonal) susceptibilities of the charm quantum number, as well as baryon–charm (BC) correlations and their ratios.

The right panel of Fig. 4 illustrates the diagonal susceptibilities $ \chi^2_\ell $, $ \chi^2_s $, and $ \chi^2_c $. At low temperatures, the quark flavor susceptibilities vanish because the chiral quark condensates have large values, resulting in large effective quark masses and suppressed fluctuations. At high temperatures, where the relevant degrees of freedom are liberated, quark masses become small, and temperature dominates the dynamics. The system behaves like a non-interacting ideal gas of massless quarks and gluons, and the susceptibilities close the corresponding SB-limits. The susceptibilities increase monotonically with temperature. The PQM model exhibits a smooth crossover before saturation, with values remaining below the SB-limit. Both $ \chi^2_\ell $ and $ \chi^2_s $ saturate at approximately 95% of the SB-limit. The light quark susceptibility $ \chi^2_\ell $ increases more rapidly than the strange and charm quark susceptibilities, consistent with the faster melting of the light chiral condensate. We also find that the pseudo-critical temperature associated with the charm quark susceptibility $ \chi^2_c $ is higher than those for the light and strange quarks, lying beyond the crossover region. It is important to note that the PQM results incorporating the Polyakov potential (Eq. (8)) improve the agreement between the PQM model and lattice QCD calculations [74].

The main objective here is to analyze how thermal effects on the charm quark mass are influenced by isospin asymmetry and to compare these results with available lattice QCD data. While the momentum dependence of the charm quark mass has been taken into account, the thermal behavior of $ \chi^2_c $ in the low-temperature region appears fixed, with noticeable evolution only at higher temperatures. Finally, the inclusion of the charm quark susceptibility in the PQM model enhances the spatial and thermal resolution of comparisons with lattice QCD simulations. Figure 5 shows the temperature dependence of baryon–charm correlations $ \chi^{\rm BC}_{nm} $ (left panel (a)), higher-order fluctuations of the charm quantum number $ \chi^C_n $ (middle panel (b)), and the corresponding ratios between correlations and fluctuations, compared with lattice QCD simulations [76] (right panel (c)), all at vanishing chemical potential $ \mu_f = 0\; \text{MeV} $. The present SU(4) PQM results (solid curves) are obtained in the presence of the Polyakov loop variables.

Figure 5.  (color online) Thermal dependence of baryon-charm correlations (left panel (a)), higher-order charm fluctuations (middle panel (b)), and their ratios between correlations and fluctuations (right panel (c)) at vanishing chemical potential, $ \mu_f=0$ MeV. The PQM results are compared with lattice QCD simulations (solid symbols) [76]. The solid arrows indicate the corresponding SB limits.

The right panel (a) of Fig. 5 shows the thermal behavior of the correlation between baryon and charm quantum numbers, denoted as $ \chi^{\rm BC}_{nm} $, where $ (n,m) $ refers to the order of the derivatives with respect to charm and baryon chemical potentials, respectively. The solid arrows in the plot indicate the expected SB- limits for each correlation. These correlations are sensitive to the interplay between baryon number and charm quantum number and are important observables for probing the liberation of charm degrees of freedom across the QCD crossover transition. At low temperatures, the correlations are small, indicating that charm degrees of freedom are still confined in the hadronic phase. As the temperature enters the crossover region, the correlations rapidly increase, signaling the partial restoration of chiral symmetry and the onset of deconfinement for charm-bearing states. At high temperatures, these correlations reach a saturation level below the SB-limit, as indicated by arrows. Their sensitivity to charm degrees of freedom is evident, as the pseudo-critical temperature shifts lower with increasing derivative order with respect to the charm quantum number. A systematic suppression in thermal behavior is also observed with higher-order derivatives. The influence of charm fluctuations becomes more pronounced with the inclusion of the charm chiral condensate and charm quantum number in the PQM model, enabling a more accurate evaluation of the quark–hadron phase structure at high temperatures. Additionally, the incorporation of charm meson states improves the spatial resolution of these comparisons, particularly in the crossover and high-temperature regions.

Different $ (n, m) $ combinations give a different response, where the lowest-order correlation $ \chi^{\rm BC}_{11} $ (black solid curve) rises most strongly and saturates close to its SB-limit. The higher-order derivatives, such as $ \chi^{\rm BC}_{13} $ (blue dashed) and $ \chi^{\rm BC}_{22} $ (red dashed), grow more slowly and saturate at lower values, while $ \chi^{\rm BC}_{31} $ (brown dash-dot) remains suppressed over the full temperature range and contributes negligibly to the SB limit. We conclude that the sensitivity of the observables to charm fluctuations increases with the inclusion of higher derivatives, making them useful probes for testing the strength of QCD interactions, especially in the transition region relevant to heavy-ion collision experiments like those at NICA and FAIR.

The middle panel (b) of Fig. 5 illustrates the higher-order fluctuations of the charm quantum number, $ \chi_n^C $, for $ n = 2,\ 4, $ and 6 in the SU(4) PQM model. In the low-temperature region, the thermal behavior of $ \chi_n^C $ shows small values, characteristic of the hadronic phase. This suppression is due to large chiral condensates, which give rise to large effective masses, thereby reducing susceptibility to fluctuations. The evolution of $ \chi_n^C $ with temperature does not exhibit a sharp phase transition but rather a smooth crossover. The temperature dependence of fluctuations is significantly influenced by the chiral structure of the hadronic states. At high temperatures, a new state of matter emerges, characterized by nearly massless quarks and gluons. As the system transitions to the deconfined phase, the charm fluctuations gradually increase and approach a saturation value slightly below the SB-limit. The higher-order fluctuations, i.e., $ \chi_6^C $ are more sensitive to the critical behavior of the system. Their response is sharper and peaks earlier than lower-order fluctuations. Notably, the transition to the partonic phase becomes more rapid as the order n of the derivative of the charm quantum number increases, indicating enhanced sensitivity of higher-order fluctuations to critical dynamics and deconfinement features.

The right panel (c) of Fig. 5 illustrates the ratios between correlations and fluctuations of baryon–charm quantum numbers, $ \chi^{\rm BC}_{nm}/\chi^{\rm BC}_{\alpha \beta} $, as functions of temperature T at vanishing chemical potential ($ \mu_f = 0\; $ MeV), where $ n, m, \alpha, \beta > 0 $ and $ n + m $, $ \alpha + \beta $ are even. The ratios $ \chi_4^C/\chi_2^C $, $ \chi_6^C/\chi_2^C $, $ \chi^{\rm BC}_{11}/\chi^{BC}_{13} $, and $ \chi^{\rm BC}_{31}/\chi^{\rm BC}_{13} $ are obtained from the SU(4) PQM model at vanishing chemical potential $ \mu_f = 0\; $MeV. The thermal behaviors of the ratios $ \chi_4^C/\chi_2^C $ and $ \chi^{\rm BC}_{11}/\chi^{\rm BC}_{13} $ are similar in strength. In the low-temperature and crossover regions, these ratios begin near unity, implying that the corresponding correlations and fluctuations are comparable in strength in the hadronic phase i.e., $ \chi_4^C \approx \chi_2^C $ and $ \chi^{\rm BC}_{11} \approx \chi^{\rm BC}_{13} $. In the hadronic phase, the PQM results agree well with lattice QCD simulations [76].

As the temperature increases and the system enters the crossover region, both ratios gradually decrease and then stabilize at higher temperatures, indicating that higher-order baryon–charm correlations become weaker in the deconfined phase. This behavior suggests that baryon– charm interactions follow a perturbative pattern at high T, with the ratios approaching values below the SB-limit. The thermal behavior of the ratio $ \chi_6^C/\chi_2^C $ also begins near unity in the crossover region, confirming that charm fluctuations $ \chi_2^C $ and $ \chi_6^C $ are nearly equal at low temperature– well-representing the contribution of charm degrees of freedom. In this region, the suppression is due to large chiral condensates, which give rise to large effective masses, leading to smaller values for higher-order susceptibilities. As temperature increases, the ratio $ \chi_6^C/\chi_2^C $ rapidly drops, reflecting a sharp change in the system near the crossover–consistent with signatures of the chiral phase transition. This decline indicates that higher-order charm fluctuations are significantly suppressed in the deconfined phase, although charm interactions remain present. The ratios $ \chi^{\rm BC}_{31}/\chi^{\rm BC}_{13} $ and $ \chi^{\rm BC}_{13}/\chi^{\rm BC}_{22} $ stabilize around 0.1 and 3.0, respectively. This consistent suppression indicates a strong reduction of higher-order baryon–charm correlations relative to lower-order ones and reflects how the system evolves from a hadronic to partonic phase. The PQM model agrees well with lattice QCD calculations in the crossover region, confirming its reliability for studying QCD thermodynamics with charm quark effects, as it accurately captures charm-related fluctuations and matches lattice QCD results at high temperatures.

The squared masses of the meson states are defined as the second derivatives of the grand canonical PQM thermodynamic free energy density (Eq. (11)) with respect to the corresponding mesonic fields at the global minimum

where the subscripts $ \sigma_{i,a(b)} $ denote the corresponding meson field of $ i- $ hadron state, $ a,b\in[0,1,\cdots,15] $. Eq. (23) is defined as a follow:

● The first term, $ m_{i,ab}^2 \big|_{\chi} $, represents the tree-level mass matrix obtained from the second derivatives of the chiral Lagrangian with respect to the meson fields. Tables 2 and 3 provide the analytical expressions for the (pseudo)scalar meson sectors, including the anomaly term in the SU(4) PQM model: The first two columns assign the meson field in terms of $ pi_i $ pseudoscalar meson field and PDG counterparts [77]. The third column gives the quark flavor content based on SU(4) symmetry. Furthermore, it gives the tree-level mass-squared derived analytically from the SU(4) PQM chiral Lagrangian in Eq. (1).

Field	PDG/MeV	Quark content	Mass squared	Analytical derivation
$ \pi^\pm=(\pi_1\pm {\rm i} \pi_2)/\sqrt{2} $, $ \pi^0=\pi_3 $	$ \pi(1300) $	$ u\bar{d}, d\bar{u} $ $ \dfrac{u\bar{u}-d\bar{d}}{\sqrt{2}} $	$ m_{\pi}^2 $	$ m_0^2 +\big(\lambda_1 +\dfrac{\lambda_2}{2} \big) \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2}\big) + \lambda_1 \big(\sigma_s^2 +\sigma_c^2\big) - \dfrac{{\cal{C}}}{2} \sigma_s \sigma_c $.
$ K^\pm=(\pi_4\pm {\rm i} \pi_5)/\sqrt{2} $, $ K^0 =(\pi_6+ {\rm i} \pi_7)/\sqrt{2} $, $ \overline{K}^0=(\pi_6 - {\rm i} \pi_7)/\sqrt{2} $	$ K(1460) $	$ u\bar{s}, s\bar{u} $ $ d\bar{s} $ $ s\bar{d} $	$ m^2_K $	$ m_0^2 + \big(\lambda_1 +\dfrac{\lambda_2}{2}\big) \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2} \big) - \lambda_2 \big(\dfrac{\sigma_u}{\sqrt{2}} + \dfrac{\sigma_d}{\sqrt{2}} \big) \sigma_s+\big(\lambda_1 + \lambda_2\big) \sigma^2_s+ \lambda_1 \sigma_c^2 - \dfrac{{\cal{C}}}{2} \sigma_s \sigma_c $.
$ \eta_\ell $	$ \eta_\ell (1408) $	$ \dfrac{u\bar{u}+d\bar{d}}{\sqrt{2}} $	$ m_{\eta_\ell}^2 $	$ m_0^2 +\big(\lambda_1 +\dfrac{\lambda_2}{2} \big) \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2}\big) + \lambda_1 \big(\sigma_s^2 +\sigma_c^2\big) + \dfrac{{\cal{C}}}{2} \sigma_s \sigma_c $.
$ \eta_s $	$ \eta_s(1475) $	$ s\bar{s} $	$ m_{\eta_s}^2 $	$ m_0^2 +\lambda_1 \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2} + \sigma_s^2 +\sigma_c^2\big) + \lambda_2 \sigma_s^2 $.
Eigenv- ($ \eta_\ell,\eta_s $)	$ \eta (1294) $	$ \dfrac{u\bar{u}-d\bar{d}}{\sqrt{2}} $ $ s\bar{s} $	$ m_{\eta}^2 $	$ \dfrac{1}{2} \Big(m_{\eta_\ell}^2 +m_{\eta_s}^2 - \Delta_{\eta_\ell,\eta_s} \Big) $.
Eigenv+ ($ \eta_\ell,\eta_s $)	$ \eta^\prime (958) $	$ \dfrac{u\bar{u}-d\bar{d}}{\sqrt{2}} $ $ s\bar{s} $	$ m_{\eta^\prime}^2 $	$ \dfrac{1}{2} \Big(m_{\eta_\ell}^2 +m_{\eta_s}^2 + \Delta_{\eta_\ell,\eta_s} \Big) $, $ \Delta_{\eta_\ell,\eta_s} =\sqrt{\big(m_{\eta_\ell}^2 - m_{\eta_s}^2 \big)^2 + 4m_{{\eta_\ell,\eta_s}}^2 } $, $ m_{\eta_\ell,\eta_s}^2 = \dfrac{{\cal{C}}}{2} \big(\dfrac{\sigma_u}{2} + \dfrac{\sigma_d}{2}\big) \sigma_c $.
$ \eta_c=(\pi_0-\sqrt{3} \pi_{15})/2 $	$ \eta_c (2983) $	$ c\bar{c} $	$ m_{\eta_c}^2 $	$ m_0^2 +\lambda_1 \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2} + \sigma_s^2 +\sigma_c^2\big) + \lambda_2 \sigma_c^2 +2 \varepsilon_c $.
$ D_0^\pm =(\pi_{11} \pm {\rm i} \pi_{12})/\sqrt{2} $, $ D_0^0 =(\pi_{9} - {\rm i} \pi_{10})/\sqrt{2} $, $ \overline{D}_0^0=(\pi_{9} + {\rm i} \pi_{10})/\sqrt{2} $	$ D^0(1865) $	$ u\bar{c},c\bar{u} $	$ m_{D}^2 $	$ m_0^2 + \big(\lambda_1 +\dfrac{\lambda_2}{2} \big) \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2}\big) + \lambda_1 (\sigma_s^2 + \sigma_c^2) + \lambda_2 \sigma_c^2 -\dfrac{\lambda_2}{\sqrt{2}} \big(\dfrac{\sigma_u}{2} + \dfrac{\sigma_d}{2}\big) \sigma_c + \varepsilon_c $.
$ D_{S,0}^\pm=(\pi_{13} \pm {\rm i} \pi_{14})/\sqrt{2} $	$ D_{S}^\pm (1968) $	$ s\bar{c}, c\bar{s} $	$ m_{D_s}^2 $	$ m_0^2 + \lambda_1 \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2} + \sigma_s^2 +\sigma_c^2\big) +\lambda_2 (\sigma_s^2 + \sigma_c^2) -\lambda_2 \sigma_s \sigma_c + \varepsilon_c $.

Table 2.  Mass expression of spin $ J^{{\cal{P}} {\cal{C}}} = 0^{-+} $ pseudoscalar meson sectors, PDG correspondence, quark content, quantum numbers, and PDG mass values from [77].

Field	PDG/MeV	Quark content	Mass squared	Analytical derivation
$ a_0^\pm=(\sigma_1\pm{\rm i}\sigma_2)/\sqrt{2} $, $ a_0^0=\sigma_3 $	$ a_0(1474) $	$ u\bar{d}, d\bar{u} $ $ \dfrac{u\bar{u}-d\bar{d}}{\sqrt{2}} $	$ m_{a_0}^2 $	$ m_0^2 +\lambda_1 \big(\sigma_u^2 +\sigma_d^2 - \sigma_u\sigma_d + \sigma_s^2 +\sigma_c^2\big) + \dfrac{3\lambda_2}{2}\big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2}\big) +\dfrac{{\cal{C}}}{2} \sigma_s \sigma_c $,
$ \kappa^\pm=(\sigma_4\pm {\rm i} \sigma_5)/\sqrt{2} $, $ \kappa^0 =(\sigma_6+ {\rm i} \sigma_7)/\sqrt{2} $, $ \overline{\kappa}^0=(\sigma_6 - {\rm i} \sigma_7)/\sqrt{2} $	$ k^*_0(1425) $	$ u\bar{s}, s\bar{u} $ $ d\bar{s} $ $ s\bar{d} $	$ m^2_\kappa $	$ m_0^2 + \big(\lambda_1 +\dfrac{\lambda_2}{2}\big) \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2} \big) +\lambda_2 \big(\dfrac{\sigma_u}{\sqrt{2}} + \dfrac{\sigma_d}{\sqrt{2}} \big) \sigma_s+\big(\lambda_1 + \lambda_2\big) \sigma^2_s+ \lambda_1 \sigma_c^2 + \dfrac{{\cal{C}}}{2} \sigma_s \sigma_c $
$ \sigma_\ell $	$ \sigma_\ell (1370) $	$ \dfrac{u\bar{u}+d\bar{d}}{\sqrt{2}} $	$ m_{\sigma_\ell}^2 $	$ m_0^2 +\lambda_1 \big(\sigma_u^2 +\sigma_d^2 + \sigma_u\sigma_d + \sigma_s^2 +\sigma_c^2\big) + \dfrac{3\lambda_2}{2}\big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2}\big) -\dfrac{{\cal{C}}}{2} \sigma_s \sigma_c $,
$ \sigma_s $	$ \sigma_s(1506) $	$ s\bar{s} $	$ m_{\sigma_s}^2 $	$ m_0^2 + \lambda_1\big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2}\big) + 3(\lambda_1+\lambda_2)\sigma s^2 + \lambda_1 \sigma_c^2, $
Eigenv- ($ \sigma_\ell,\sigma_s $)	$ f_0(1370) $	$ \dfrac{u\bar{u}+d\bar{d}}{\sqrt{2}} $ $ s\bar{s} $	$ m_{f_0}^2 $	$ \dfrac{1}{2} \Big(m_{\sigma_\ell}^2 +m_{\sigma_s}^2 - \Delta_{\sigma_\ell,\sigma_s} \Big) $.
Eigenv+ ($ \sigma_\ell,\sigma_s $)	$ \sigma (1200) $	$ \dfrac{u\bar{u}+d\bar{d}}{\sqrt{2}} $ $ s\bar{s} $	$ m_{\sigma}^2 $	$ \dfrac{1}{2} \Big(m_{\sigma_\ell}^2 +m_{\sigma_s}^2+ \Delta_{\sigma_\ell,\sigma_s} \Big) $,V$ \Delta_{\sigma_\ell,\sigma_s} =\sqrt{\big(m_{\sigma_\ell}^2 - m_{\sigma_s}^2 \big)^2 + 4m_{{\sigma_\ell,\sigma_s}}^2 } $, $ m_{\sigma_\ell,\sigma_s}^2 = \lambda_1 \big(\dfrac{\sigma_u}{2} + \dfrac{\sigma_d}{2}\big) \sigma_s - \dfrac{{\cal{C}}}{2} \big(\dfrac{\sigma_u}{2} + \dfrac{\sigma_d}{2}\big)\sigma_c $..
$ \chi_{c_0}=(\sigma_3-\sqrt{3} \pi_{15})/2 $	$ \chi_{c_0} (3415) $	$ c\bar{c} $	$ m_{\chi_{c_0}}^2 $	$ m_0^2 +\lambda_1 \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2} + \sigma_s^2 \big) + 3 (\lambda_1+ \lambda_2) \sigma_c^2 +2 \varepsilon_c $,
$ D_0^\pm =(\sigma_{11} \pm {\rm i} \sigma_{12})/\sqrt{2} $, $ D_0^0 =(\sigma_{9} - {\rm i} \sigma_{10})/\sqrt{2} $, $ \overline{D}_0^0=(\sigma_{9} + {\rm i} \sigma_{10})/\sqrt{2} $	$ D_0^*(2350) $	$ u\bar{c},c\bar{u} $	$ m_{D_0^*}^2 $	$ m_0^2 + \big(\lambda_1 +\dfrac{\lambda_2}{2} \big) \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2}\big) + \lambda_1 (\sigma_s^2 + \sigma_c^2) + \lambda_2 \sigma_c^2 +\dfrac{\lambda_2}{\sqrt{2}} \big(\dfrac{\sigma_u}{2} + \dfrac{\sigma_d}{2}\big) \sigma_c + \varepsilon_c $.
$ D_{S,0}^\pm=(\sigma_{13} \pm {\rm i} \sigma_{14})/\sqrt{2} $	$ D_{s0}*(2318) $	$ s\bar{c}, c\bar{s} $	$ m_{D_{s0}*}^2 $	$ m_0^2 + \lambda_1 \big(\dfrac{\sigma_u^2}{2} + \dfrac{\sigma_d^2}{2} + \sigma_s^2 +\sigma_c^2\big) +\lambda_2 (\sigma_s^2 + \sigma_c^2) + \lambda_2 \sigma_s \sigma_c + \varepsilon_c $.

Table 3.  Same as Table 2 but for spin $ J^{{\cal{P}} {\cal{C}}} = 0^{++} $ scalar meson sectors.

In Table 2, the squared mass of the K meson state corresponds to the ($ 4,4 $) element of the pseudoscalar mass matrix and satisfies the Gell-Mann–Okubo mass relation [78–80]. The η and $ \eta^\prime $ masses result from the diagonalization of the isoscalar pseudoscalar mass matrix, which incorporates off-diagonal flavor mixing and anomaly terms [21]. The charm isosinglet $ \eta_c $ includes $ \sigma_c $ contributions and is further shifted due to charm-sector symmetry breaking. The D resonances include charm-specific symmetry-breaking terms ($ \epsilon_c $). A similar approach applies to Table 3 but for scalar meson states. The κ meson state corresponds to the ($ 4,4 $) element of the scalar mass matrix and also satisfies the Gell-Mann–Okubo mass relation [78–80]. The isoscalar states $ f_0 $ and σ in the scalar multiplet are obtained via the diagonalization of the isoscalar scalar mass matrix, which includes off-diagonal flavor mixing and anomaly terms, as presented in Table 3 [21].

● The second term $ m_{i,ab}^2 \big|_{\text{the}} $ arises from the fermionic (partonic) contribution to the grand thermodynamic free energy, which describes the thermal and dense QCD-medium effects on meson sectors due to the quark–antiquark contribution, as given in Eq. (13) within the SU(4) PQM model [21, 22].

where $ {\cal{N}}_c = 2\,N_c $ color flavors. This expression includes the squared quark masses with respect to meson fields from the SU(4) PQM model (Eq. (24)). The diagonalizations of the quark-mass matrix are listed in Table 4. The quark mass derivative with respect to the meson fields $ \sigma_{i,a} $ is defined as $ m^2_{f,a} \equiv \partial m^2_f/\partial \sigma_{i,a} $, and the second derivative with respect to $ \sigma_{i,a} $ and $ \sigma_{i,b} $ is given by $ m^2_{f,{ab}} \equiv \partial^2 m^2_f/\partial \sigma_{i,a} \partial \sigma_{i,b} $. More details can be found in Refs. [22, 81]. Table 4 lists the squared quark masses in the SU(4) PQM model with respect to the meson fields: the first four columns correspond to the light (u and d) quarks, the fifth and sixth columns to the strange (s) quark, and the seventh and eighth columns to the charm (c) quark. The terms $ E_{q,f}(T,\mu) $ are symmetric under $ \mu \rightarrow -\mu $, i.e., $ E_{q,f}(T,\mu) = E_{q,f}(T,-\mu) $, and thus ensure that the thermodynamic contributions from quark–antiquark pairs remain even functions of the chemical potential. This symmetry plays an essential role in conserving CP invariance and simplifies the evaluation of thermodynamic integrals at vanishing and finite μ:

Field	$ \dfrac{m_{u,\alpha}^2 m_{u,\beta}^2}{g^4} $	$ \dfrac{m_{u,\alpha\beta}^2 }{g^2} $	$ \dfrac{m_{d,\alpha}^2 m_{d,\beta}^2}{g^4} $	$ \dfrac{m_{d,\alpha\beta}^2}{g^2} $	$ \dfrac{m_{s,\alpha}^2 m_{s,\beta}^2}{g^4} $	$ \dfrac{m_{s,\alpha\beta}^2}{g^2} $	$ \dfrac{m_{c,\alpha}^2 m_{c,\beta}^2 }{g^4} $	$ \dfrac{m_{c,\alpha\beta}^2}{g^2} $
$ \sigma_{0} \sigma_{0} $	$ \dfrac{1}{8} \sigma_u^2 $	$ \dfrac{1}{4} $	$ \dfrac{1}{8} \sigma_d^2 $	$ \dfrac{1}{4} $	$ \dfrac{1}{4} \sigma_s^2 $	$ \dfrac{1}{4} $	$ \dfrac{1}{4} \sigma_c^2 $	$ \dfrac{1}{4} $
$ \sigma_0 \sigma_3 $	$ \dfrac{\sigma_u^2}{4\sqrt{2}} $	$ \dfrac{1}{2\sqrt{2}} $	$ -\dfrac{\sigma_d^2}{4\sqrt{2}} $	$ - \dfrac{1}{2\sqrt{2}} $	0	0	0	0
$ \sigma_{0} \sigma_{8} $	$ \dfrac{1}{4\sqrt{6}} \sigma_u^2 $	$ \dfrac{1}{2\sqrt{6}} $	$ \dfrac{1}{4\sqrt{6}} \sigma_d^2 $	$ \dfrac{1}{2\sqrt{6}} $	$ -\dfrac{\sqrt{2}}{3} \sigma_s^2 $	$ \dfrac{1}{\sqrt{6}} $	0	0
$ \sigma_0 \sigma_{15} $	$ \dfrac{1}{8\sqrt{3}}\sigma_u^2 $	$ \dfrac{1}{4\sqrt{3}} $	$ \dfrac{1}{8\sqrt{3}}\sigma_d^2 $	$ \dfrac{1}{4\sqrt{3}} $	$ -\dfrac{1}{4\sqrt{3}}\sigma_s^2 $	$ \dfrac{1}{4\sqrt{3}} $	$ \dfrac{\sqrt{3}}{4} \sigma_c^2 $	$ -\dfrac{\sqrt{3}}{4} $
$ \sigma_{3} \sigma_{3} $	$ \dfrac{1}{4} \sigma_u^2 $	$ \dfrac{1}{2} $	$ -\dfrac{1}{4} \sigma_d^2 $	$ -\dfrac{1}{2} $	0	0	0	0
$ \sigma_3 \sigma_8 $	$ \dfrac{1}{4\sqrt{3}} \sigma_u^2 $	$ \dfrac{1}{2\sqrt{3}} $	$ -\dfrac{1}{4\sqrt{3}} \sigma_u^2 $	$ -\dfrac{1}{2\sqrt{3}} $	0	0	0	0
$ \sigma_3 \sigma_{15} $	$ \dfrac{1}{4\sqrt{6}} \sigma_u^2 $	$ \dfrac{1}{2\sqrt{6}} $	$ -\dfrac{1}{4\sqrt{6}} \sigma_u^2 $	$ -\dfrac{1}{2\sqrt{6}} $	0	0	0	0
$ \sigma_{8} \sigma_{8} $	$ \dfrac{1}{12} \sigma_u^2 $	$ \dfrac{1}{6} $	$ \dfrac{1}{12} \sigma_d^2 $	$ \dfrac{1}{6} $	$ \dfrac{2}{3} \sigma_s^2 $	$ \dfrac{2}{3} $	0	0
$ \sigma_8 \sigma_{15} $	$ \dfrac{1}{12\sqrt{2}} \sigma_u^2 $	$ \dfrac{\sqrt{1}}{6\sqrt{2}} $	$ \dfrac{1}{12\sqrt{2}} \sigma_d^2 $	$ \dfrac{\sqrt{1}}{6\sqrt{2}} $	$ \dfrac{1}{3\sqrt{2}} \sigma_s^2 $	$ \dfrac{\sqrt{1}}{3\sqrt{2}} $	0	0
$ \sigma_{15} \sigma_{15} $	$ \dfrac{1}{24} \sigma_u^2 $	$ \dfrac{1}{12} $	$ \dfrac{1}{24} \sigma_d^2 $	$ \dfrac{1}{12} $	$ \dfrac{1}{12} \sigma_s^2 $	$ \dfrac{1}{12} $	$ \dfrac{3}{4} \sigma_c^2 $	$ \dfrac{3}{4} $
$ \sigma_{1} \sigma_{1} $	$ \sigma_u^2 $	1	$ \sigma_d^2 $	1	0	0	0	0
$ \sigma_{4} \sigma_{4} $	0	$ \sigma_u \dfrac{\sigma_u + \sqrt{2} \sigma_s}{\sigma_u^2 - 2 \sigma_s^2} $	0	$ \sigma_d \dfrac{\sigma_d + \sqrt{2} \sigma_s}{\sigma_d^2 - 2 \sigma_s^2} $	0	$ \sigma_s \dfrac{2\sigma_s + \sqrt{2} \sigma_{\ell}}{2 \sigma_s^2 - \sigma_{\ell}^2} $	0	0
$ \sigma_{6} \sigma_{6} $	0	0	0	0	0	$ \sigma_s \dfrac{\sigma_s + \sqrt{3} \sigma_{c}}{\sigma_s^2 -3\sigma_{c}^2} $	0	$ \sigma_c \dfrac{\sqrt{3} \sigma_s + 3 \sigma_{c}}{3\sigma_c^2 -\sigma_{s}^2} $
$ \pi_{0} \pi_{0} $	0	$ \dfrac{1}{4} $	0	$ \dfrac{1}{4} $	$ \dfrac{1}{4} \sigma_s^2 $	$ \dfrac{1}{4} $	$ \dfrac{1}{4} \sigma_c^2 $	$ \dfrac{1}{4} $
$ \pi_0 \pi_3 $	0	$ \dfrac{1}{2\sqrt{2}} $	0	$ - \dfrac{1}{2\sqrt{2}} $	0	0	0	0
$ \pi_{0} \pi_{8} $	0	$ \dfrac{1}{2\sqrt{6}} $	0	$ \dfrac{1}{2\sqrt{6}} $	$ -\dfrac{\sqrt{2}}{3} \sigma_s^2 $	$ \dfrac{1}{\sqrt{6}} $	0	0
$ \pi_0 \pi_{15} $	0	$ \dfrac{1}{4\sqrt{3}} $	0	$ \dfrac{1}{4\sqrt{3}} $	$ -\dfrac{1}{4\sqrt{3}}\sigma_s^2 $	$ \dfrac{1}{4\sqrt{3}} $	$ \dfrac{\sqrt{3}}{4} \sigma_c^2 $	$ -\dfrac{\sqrt{3}}{4} $
$ \pi_{3}\pi_{3} $	0	$ \dfrac{1}{2} $	0	$ -\dfrac{1}{2} $	0	0	0	0
$ \pi_3 \pi_8 $	0	$ \dfrac{1}{2\sqrt{3}} $	0	$ -\dfrac{1}{2\sqrt{3}} $	0	0	0	0
$ \pi_3 \pi_{15} $	0	$ \dfrac{1}{2\sqrt{6}} $	0	$ -\dfrac{1}{2\sqrt{6}} $	0	0	0	0
$ \pi_{8} \pi_{8} $	0	$ \dfrac{1}{6} $	0	$ \dfrac{1}{6} $	0	$ \dfrac{2}{3} $	0	0
$ \pi_8 \pi_{15} $	0	$ \dfrac{\sqrt{1}}{6\sqrt{2}} $	0	$ \dfrac{\sqrt{1}}{6\sqrt{2}} $	0	$ \dfrac{\sqrt{1}}{3\sqrt{2}} $	0	0
$ \pi_{15} \pi_{15} $	0	$ \dfrac{1}{12} $	0	$ \dfrac{1}{12} $	0	$ \dfrac{1}{12} $	0	$ \dfrac{3}{4} $
$ \pi_{1} \pi_{1} $	0	1	0	1	0	0	0	0
$ \pi_{4} \pi_{4} $	0	$ \sigma_u \dfrac{\sigma_u - \sqrt{2} \sigma_s}{\sigma_u^2 - 2 \sigma_s^2} $	0	$ \sigma_d \dfrac{\sigma_d - \sqrt{2} \sigma_s}{\sigma_d^2 - 2 \sigma_s^2} $	0	$ \sigma_s \dfrac{2\sigma_s - \sqrt{2} \sigma_{\ell}}{2 \sigma_s^2 - \sigma_{\ell}^2} $	0	0
$ \pi_{6} \pi_{6} $	0	0	0	0	0	$ \sigma_s \dfrac{\sigma_s - \sqrt{3} \sigma_{c}}{\sigma_s^2 -3\sigma_{c}^2} $	0	$ \sigma_c \dfrac{\sqrt{3} \sigma_s + 3 \sigma_{c}}{3\sigma_c^2 -\sigma_{s}^2} $

Table 4.  Squared–quark mass second derivatives with respect to the relevant meson fields within SU(4) PQM model.

are implemented as in Ref. [81]. Furthermore, for (anti-) quark, $ {\cal{B}}_{q_f (\bar{q}_f)} =3 {\cal{N}}_{q_f (\bar{q}_f)}^2 - {\cal{C}}_{q_f (\bar{q}_f)} $, where

are defined as in Ref. [81].

● The third term $ m_{i,ab}^2 \big|_{\text{vac}} $ defines the fermionic vacuum correction [82, 83]. Using the three-dimensional– cutoff Λ in the fermionic vacuum term in Eq. (13), one obtains

The present calculations of meson state squared masses account for both the isospin asymmetry of light quarks ($ m_u \neq m_d $) and the inclusion of heavy quark masses ($ m_s $ and $ m_c $). Table 5 provides a comprehensive comparison of the (pseudo)scalar meson masses in the SU(3) (third column) and SU(4) (sixth column) formulations of the PQM model, both incorporating isospin asymmetry. Our results are compared with those of other effective thermal QCD models, including the extended Linear sigma model (eLSM) [30, 31] and PNJL model [84], as well as the Particle Data Group (PDG) values [77]. The PDG column shows the experimental meson masses, with ground states listed first and their radial excitations added where available. Radial excitations are simply heavier partners of the ground states, carrying the same quantum numbers but reflecting a higher internal quark–antiquark mode. The results exhibit strong consistency with experimental data and with predictions from other QCD-based effective models. In particular, the agreement between the SU(4) PQM model and eLSM is notably fine, especially in both the scalar and pseudoscalar meson sectors. The extended SU(4) PQM model not only reproduces the ground-state spectrum but also matches well with several known excitations, such as $ \pi (1300), \eta(1405/1475) $, and $ f_0(1370) $. This agreement suggests that the model is able to capture important aspects of chiral symmetry restoration and confinement dynamics across the light, strange, and charm sectors.

Meson	PDG [77] (ground; excitation)	SU(3) PQM	eLSM [30, 31]	PNJL [84]	SU(4) PQM
Pseudoscalar Sector ($ 0^{-+} $)
π	$ \pi^0(134.98 \pm 0.06) $; $ \pi(1300 \pm 100) $	136.8	$ 141.0\pm5.8 $	126	1291
K	$ K^0(497.61 \pm 0.01) $; $ K(1460) $	496.3	$ 485.6\pm3.0 $	490	1228
$ \eta_\ell $	–	–	–	–	1342
$ \eta_s $	–	–	–	–	1469
η	$ \eta(547.86 \pm 0.02) $; $ \eta(1294 \pm 4) $	553	$ 509.4\pm3.0 $	505	1276
$ \eta^\prime $	$ \eta^\prime(957.78 \pm 0.06) $; $ \eta^\prime(1475 \pm 4) $	635	$ 962.5\pm 5.6 $	949	1527
$ \eta_c $	$ \eta_c(2983.9 \pm 0.4) $	–	$ 2673\pm 118 $	–	2682
D	$ D^0(1864.84 \pm 0.05) $	–	$ 1981\pm73 $	–	1985.73
$ D_s $	$ D_s^\pm(1968.35 \pm 0.07) $	–	$ 2004\pm73 $	–	2035.02
Scalar Sector ($ 0^{++} $)
$ a_0 $	$ a_0(980 \pm 20) $; $ a_0(1474 \pm 19) $	849.8	$ 1636\pm1 $	837	1381.58
κ	$ k^*_0(682 \pm 29) $; $ k^*_0(1425 \pm 50) $	1331.6	$ 1450\pm 1 $	1013	1581.45
$ \sigma_\ell $	–	–	–	–	1458.01
$ \sigma_s $	–	–	–	–	1671.23
σ	400–1200; $ \sigma(1704 \pm 12) $	824.0	$ 975.1\pm 6.4 $	700	1645.96
$ f_0 $	1200–1500	1339.0	$ 1186 \pm 6 $	1169	1241.92
$ \chi_{c_0}(1P) $	$ \chi_{c_0}(3414.71 \pm 0.3) $	–	$ 3144\pm128 $	–	3541.21
$ D_0^* $	$ D_0^*(2343.00 \pm 10.0) $	–	$ 2414\pm77 $	–	2451.11
$ D_{s0}^* $	$ D_{s_0}^* (2317.8 \pm 0.5) $	–	$ 2467\pm76 $	–	2507.47

Table 5.  Comparison between (pseudo)scalar meson states within PQM model.

Table 5 aims to estimate meson state masses with high precision, not only for vacuum phenomenology but also for applications in thermal and dense QCD media. The available fits in the SU(3) and SU(4) PQM models are obtained from the expressions for the squared tree-level masses listed in Tables 2 and 3, corresponding to the pseudoscalar and scalar sectors, respectively. The estimated fits of nonet meson states in the SU(3) framework assume a vanishing charm condensate $ \sigma_c $. The SU(3) PQM model accurately predicts the π and K meson masses, showing excellent agreement with PDG values (within 1%–2%), thereby confirming the effectiveness of the chiral Lagrangian in describing the Nambu–Goldstone bosons. However, it significantly underestimates the η and $ \eta^\prime $ masses, revealing limitations in capturing the $ U(1)_A $ anomaly within the SU(3) framework. Moreover, the states $ \eta_\ell,\; \eta_s,\; \sigma_\ell $, and $ \sigma_s $ are model basis states (light- and strange-quark isosinglets). The predicted masses of $ \eta_\ell(1342) $ and $ \eta_s(1469) $ correspond closely to the experimental $ \eta(1408) $ and $ \eta(1475) $ states, respectively, while the $ \sigma_\ell (1458.01) $ and $ \sigma_s (1671.23) $ states correspond with the experimental $ \sigma (1370) $ and $ \sigma(1506) $ states, respectively.

In the SU(4) PQM model, the focus shifts to the neutral multiplet masses of nonstrange–charmed meson states. This approach employs a low-energy chiral model in which the input parameters directly determine the calculated meson masses. The extended flavor symmetry enables exploration of charm-sector phenomenology while maintaining consistency with the light-quark sector. The choice of fitting parameters plays an essential role in achieving agreement with PDG values and results from other low-energy effective models. In the present work, the parameters are fixed at $ \sigma = 800 $ MeV, as specified in Ref. [21]. In the presence of chiral symmetry breaking, isospin asymmetry among light quarks, the anomaly term, and the Polyakov loop corrections, Fig. 6 presents the masses of nonet meson sectors: pseudoscalar (negative parity) and scalar (positive parity) meson states (left and middle panels), along with open and hidden charmed mesons in the thermal QCD medium (right panel).

Figure 6.  (color online) In-medium meson masses for the nonet pseudoscalar (left panel), scalar (middle panel), and open/hidden charmed meson sectors (right panel).

Figure 6 illustrates the thermal evolution of the meson mass spectrum for the nonet pseudoscalar mesons (π, K, $ \eta_\ell $, $ \eta_s $, η, $ \eta^\prime $) (left panel), nonet scalar mesons ($ a_0 $, κ, $ \sigma_\ell $, $ \sigma_s $, $ f_0 $, σ) (middle panel), and open/hidden charmed mesons ($ D^0 $, $ D^\pm $, $ D_0^* $, $ D_{s0}^* $, $ \eta_c $, $ \chi_{c0} $) (right panel), all given at vanishing chemical potential $ \mu_f $. The temperature dependence of the meson mass spectrum can be categorized into three regions: bosonic (mesonic), phase transition, and fermionic (partonic) regions. In the bosonic region, meson masses decrease gradually with increasing temperature, reflecting the dominance of thermal suppression of mesonic degrees of freedom. The phase transition region is characterized by a smooth but rapid decline in mass values, corresponding to the chiral crossover. At high temperatures, in the fermionic region, quark–antiquark contributions dominate, leading to a moderate increase in the meson masses due to the restoration of chiral symmetry.

At low temperatures ($ T < T_{\chi} $), bosonic contributions dominate. The tree-level meson masses presented in Tables 2 and 3 effectively restore the mass gap characteristic of the chiral phase structure in meson states. As previously discussed, the system remains in the confined phase, characterized by large chiral condensates $ \sigma_f $ and small values of the Polyakov loop order parameters. This regime allows for the estimation of vacuum mass values for each meson sector at vanishing temperature ($ T = 0 $ MeV). In this region, meson masses are expressed as functions of pure mesonic condensates $ \sigma_f $. The thermal behavior of meson masses remains nearly constant, relative to their corresponding vacuum values, until the temperature approaches the pseudo-critical point $ T \sim T_{\chi} $. It is evident that in this low-temperature regime, the masses of all meson states are effectively independent of temperature.

$ U(1)_A $ symmetry breaking appears in the tree-level meson masses through the anomaly term $ {\cal{C}} $, which is particularly significant for nonet meson states (left and middle panels) due to strange and charm quark contributions. In the SU(4) PQM model, the phase structure of the charm chiral condensate $ \sigma_c $ plays an essential role in determining the chiral behavior of meson masses. When $ \sigma_c $ vanishes, the SU(3) structure is effectively restored. In this limit, the $ U(1)_A $ anomaly effect disappears in the charm–strange meson sectors, and the term $ {\cal{L}}_{{\rm{emass}}} $ becomes dominant. This term is directly related to the explicit quark mass contribution. For simplicity, we assume $ \varepsilon_u = \varepsilon_d = \varepsilon_s = 0 $ and $ \varepsilon_c \sim m_c $.

In the phase transition region, beyond the chiral temperature $ T \sim T_{\chi} $, meson masses gradually decrease with increasing temperature. This region exhibits a smooth decline, characteristic of a slow crossover transition. Beyond the chiral temperature ($ T \sim T_{\chi} $), meson masses gradually decrease with rising temperature. This behavior reflects the nature of a smooth crossover transition. Notably, the location of the chiral critical endpoint (CEP) is not universal across all meson states; rather, it depends on the specific PQM order parameters and input conditions of the QCD medium.

At high temperatures ($ T > T_{\chi} $), fermionic (partonic) thermal contributions significantly influence the thermal evolution of meson states. As the temperature increases, these contributions become dominant and lead to the degeneration of meson masses. In contrast, their impact remains negligible at low temperatures. The observed degeneracy in meson masses arises from thermal fluctuations acting on the chiral condensates $ \sigma_f $, particularly the strange $ \sigma_s $ and charm $ \sigma_c $ condensates, and reflects the process of chiral symmetry restoration [21]. These thermal fluctuations melt the light quark condensates ($ \sigma_u $, $ \sigma_d $) more rapidly than the strange and charm condensates. Consequently, nonet pseudoscalar mesons (left panel) enter the deconfined (partonic) phase earlier than nonet scalar and charmed mesons. At very high temperatures ($ T \sim 1\; {\rm{GeV}} $), the thermal evolution of the charm condensate completes the transition for strange–charm meson states. In this regime, charmed meson masses remain relatively stable until the temperature becomes extremely high. Accordingly, charmed meson states are expected to undergo phase transition and mass degeneracy as the density of the system increases. This transition enhances thermal energy to overcome the gap associated with the Fermi surface of the meson states [22, 40]. The latter differs from one meson state to another.

In conclusion, we extended the SU(4) PQM model by incorporating isospin asymmetry to explore the chiral phase structure of QCD matter in a thermal medium. This study examined how the system undergoes phase transitions in the presence of isospin asymmetry between light quark condensates and the charm condensate. The results reveal key features of chiral symmetry restoration, deconfinement transitions, quark and charm fluctuations and correlations, thermodynamic observables, and the thermal evolution of meson masses within the SU(4) PQM framework.

The introduction of isospin asymmetry leads to several modifications in the chiral phase structure and meson dispersion relations. At finite isospin chemical potential ($ \mu_I $), the u- and d-quark condensates ($ \sigma_u $, $ \sigma_d $) split due to isospin symmetry breaking. In the SU(4) framework, all diagonal sigma-fields ($ \sigma_0 $, $ \sigma_3 $, $ \sigma_8 $, $ \sigma_{15} $) and symmetry-breaking parameters ($ h_0 $, $ h_3 $, $ h_8 $, $ h_{15} $) become nonzero, explicitly breaking SU(2) isospin symmetry through the presence of $ \sigma_3 $ and $ h_3 $. As a result, the pseudo-critical temperatures for $ \sigma_u $ and $ \sigma_d $ decrease with increasing $ \mu_I $.

The inclusion of the charm chiral condensate ($ \sigma_c $) significantly enhances the thermal resolution of the PQM model, particularly at high temperatures. While light and strange quark condensates melt rapidly near the pseudo-critical temperature, the charm condensate remains largely stable up to much higher temperatures. This enables better characterization of the crossover and high-temperature phases. Incorporating the isospin asymmetry allows the $ (T/t_chi, \; \mu_I/m_\pi) $ plane to be traced, where the pseudo-critical temperature $ T_\chi $ decreases with increasing normalized isospin chemical potential. The model agrees well with lattice QCD simulations for $ \mu_I/m_\pi < 1.0 $, and the inclusion of $ \sigma_c $ enables the exploration of high-temperature and high-density regimes.

Thermodynamic observables-such as pressure, trace anomaly, and entropy density–calculated from the SU(4) PQM model show excellent agreement with lattice QCD results. All quantities exhibit a smooth crossover and remain below the SB limit at high temperatures. These results are critical for modeling the EoS for dense QCD matter relevant to compact astrophysical objects. While $ \sigma_c $ may play a minor role in early-stage deconfinement, it could influence the EoS in the extreme interior of hybrid stars. The analysis of diagonal susceptibilities $ \chi_f^2 $ for $ f = (\ell, s, c) $ provides insight into the thermal behavior of the corresponding chiral condensates. Light and strange quark susceptibilities rise sharply near $ T_\chi $, while the charm susceptibility increases gradually, indicating delayed thermal response in the charm sector.

Meson masses were obtained from the second derivative of the SU(4) PQM grand potential with respect to the mesonic fields. This includes tree-level, thermal, and vacuum contributions, allowing for a consistent description of the meson mass spectrum. Light pseudoscalar mesons such as π, K, and η show rapid mass suppression with increasing temperature, while scalar mesons such as $ a_0 $, κ, and $ f_0 $ exhibit greater thermal stability. The inclusion of $ \sigma_c $ enables analysis of the thermal behavior of both open-charm mesons (e.g., D, $ D_s $) and hidden charmonia (e.g., $ \eta_c $, $ \chi_{c0} $), which remain thermally stable over a wide temperature range.

In dense astrophysical systems, such as neutron stars and hybrid stars, where both baryon density and temperature are extreme, charm mesons may contribute to the EoS. The slow melting of $ \sigma_c $ and the thermal robustness of hidden charmed mesons suggest potential relevance for transport and stability properties in such studies. Overall, the SU(4) PQM model provides a good framework for studying the QCD phase diagram, including heavy-flavor dynamics. It offers predictive power for both heavy-ion experiments at FAIR and NICA and for modeling extreme conditions in compact stars.

The author gratefully acknowledges the anonymous reviewer for their invaluable feedback and constructive suggestions, which have significantly improved the quality of this work.

The generator operator $ \hat{T}_a = \hat{\lambda}_a / 2 $ in U(4) is obtained from the generalized Gell-Mann matrices $ \hat{\lambda}_a $ [52], where the indices $ a = 0, \ldots, 15 $ [49]. The matrices are defined as a follow:

These generators satisfy the U(4) algebra with the commutation and anticommutation relations

where $ f_{abc} $ and $ d_{abc} $ are the antisymmetric and symmetric structure constants of SU(4), respectively. The symmetric constants are explicitly defined by

In the framework of PCAC, the decay constant $ f_a $ is related to the symmetric constants and sigma fields as

For charged and neutral pion meson ($ f_{\pi^{\pm}} =f_1,\;f_{\pi^0} =f_3 $) and kaon meson ($ f_{K^{\pm}} =f_4,\;f_{K^0} =f_6 $), the decay constants take the following form:

Hence, the isospin field $ \bar{\sigma}_3 $ can be expressed as

The scalar meson generator matrix $ T_a \sigma_a $ takes the following $ 4 \times 4 $ matrix form:

Similarly, the pseudoscalar meson generator matrix $ T_a \pi_a $ becomes: