# Thermodynamic properties and phase diagram of quark matter within non-
extensive Polyakov chiral SU (3) quark mean field model

Dhananjay Singh<EMAIL_ADDRESS>Department of Physics, Dr. B R Ambedkar
National Institute of Technology Jalandhar, Jalandhar – 144008, Punjab, India
Arvind Kumar<EMAIL_ADDRESS>Department of Physics, Dr. B R Ambedkar
National Institute of Technology Jalandhar, Jalandhar – 144008, Punjab, India

###### Abstract

In the present work, we apply Tsallis non-extensive statistics to study the
thermodynamic properties and phase diagram of quark matter in the Polyakov
chiral SU(3) quark mean field model. Within this model, the properties of the
quark matter are modified through the scalar fields
$\sigma,\zeta,\delta,\chi$, the vector fields $\omega,\rho$, $\phi$, and the
Polyakov fields $\Phi$ and $\bar{\Phi}$ at finite temperature and chemical
potential. Non-extensive effects have been introduced through a dimensionless
parameter $q$ and the results are compared to the extensive case
($q\rightarrow 1$). In the non-extensive case, the exponential in the Fermi-
Dirac (FD) function is modified to a $q$-exponential form. The influence of
$q$ parameter on the thermodynamic properties: pressure, energy, and entropy
density as well as trace anomaly is investigated. The speed of sound and
specific heat with non-extensive effects is also studied. Furthermore, the
effect of non-extensivity on the deconfinement phase transition as well as the
chiral phase transition of $u,d,$ and $s$ quarks is explored. We found that
the critical end point (CEP), which defines the point in the $(T-\mu)$ phase
diagram where the order of the phase transition changes, shifts to a lower
value of temperature, $T_{CEP}$, and a higher value of chemical potential,
$\mu_{CEP}$, as the non-extensivity is increased, i.e., $q>$1.

##  I Introduction

Quantum Chromodynamics (QCD) is a theory of strong interactions among the
fundamental constituents of matter: quarks and gluons. In the low energy
regime ($10^{-2}$ GeV), perturbative approaches to solve QCD become obsolete
and we have to rely on non-perturbative methods. This led to the development
of two theoretical approaches to QCD, namely lattice quantum chromodynamics
(LQCD) and effective models of QCD. In LQCD, numerical techniques are applied
on a discretized lattice of space and time. The discretization of space and
time allows the equations to become computationally tractable to calculate the
properties of hadrons and other QCD systems. Effective field models (EFMs)
based on effective field theory were also developed to investigate the low-
energy regime of QCD. EFM uses a simpler set of degrees of freedom than those
of the underlying theory. Studying thermodynamic properties and phase
transition of quark matter at finite temperatures has always been a topic of
great interest in high-energy physics. They play an important role in
exploring physics shortly after the Big Bang big . They are also relevant to
studying the properties and structure of compact stars at high baryonic
density and low temperature dex2014 ; baym ; li2019 ; hinderer ; weber . At
low baryonic density and high temperature, they are relevant to explain Quark
Gluon Plasma (QGP) formation in relativistic heavy-ion collision marty .
Exploring the thermodynamic properties of dense quark matter is essential to
understanding the nature of the phase transition in the ($T-\mu_{B}$) plane.
The phase diagram of QCD contains information about the equilibrium phases of
QCD and the physics of phase transition.

Many theoretical models have been used to study the properties of hot and
dense quark matter. These include Quark-Meson Coupling (QMC) model tsushima ;
saito , Polyakov QMC (PQMC) model schaefer ; herbst , Nambu-Jona-Lasinio (NJL)
model menezes ; stiele , Polyakov NJL (PNJL) model fukushima ; ratii , Linear
Sigma Model (LSM) abu , Polyakov LSM (PLSM) gatto ; schaefer2009 ; kovacs ;
kovacs2022 , Chiral SU(3) Quark-Mean-Field (CQMF) model wang ; wang2001 , and
Polyakov CQMF (PCQMF) model manisha ; nisha . However, the statistical
approach used in all these studies is often the standard Boltzmann-Gibbs (BG)
statistics. The validation of BG statistics is limited to systems that satisfy
the following conditions: short-range type interactions, homogeneous and
infinite heat bath, weak or no correlations between elements, and smooth
boundary conditions. Violation of any or more of these conditions questions
the validation of BG statistics.

Thus, Tsallis proposed non-extensive statistics which depict a more realistic
environment in heavy ion collision experiments. This formalism is an extension
of the BG statistics tsallis . In this framework the usual exponential
function is replaced by a $q-$exponential function zhao2020 ; rozy2016 ,

$exp_{q}(x)=[1+(1-q)x]^{\frac{1}{1-q}}$ (1)

where $q$ is a dimensionless parameter called the non-extensive parameter and
it contains all possible effects violating the standard BG statistics.
Correspondingly, the logarithmic function also gets modified to

$ln_{q}(x)=\frac{x^{1-q}-1}{1-q}.$ (2)

Both these functions return back to their usual forms when $q\rightarrow 1$,
i.e., $exp_{q}(x)\rightarrow exp(x)$, $ln_{q}(x)\rightarrow ln(x)$ and Tsallis
statistics returns to BG statistics. As the $q$ parameter enters into the
relevant formulae of the specific dynamical model it enables the
phenomenological examination of the stability of the model under consideration
against potential departures from the BG approach. It must be emphasized here
that Tsallis statistics is not an alternative but a generalization of the BG
statistics.

The motivation for this research comes from recent non-extensive statistics
studies which showed the violation of BG statistics in high energy physics.
There is considerable agreement between experimental data and theoretical
models based on non-extensive effects. As shown in bhat ; asmaa , BG
statistics failed to reproduce the freeze-out parameters. Tsallis statistics
has also been proven to accurately produce the transverse momentum
distributions which follow a power law distribution rather than an exponential
distribution as in the BG statistics wilk ; li ; marques ; de ; ryb . Ref.
combe  shows that Tsallis statistics can be used to explain the transport
behaviour of cold atoms in optical lattices. Tsallis statistics also describes
the behaviour of a system where ergodicity breaks down tirnakli ; cirto .
Experimental data produced by PHENIX phenix , STAR star  as well as RHIC
collaborations by ATLAS atlas , ALICE alice , and collaborations at the LHC by
the CMS cms  matches very well within the framework of Tsallis statistics.
Several authors have summarized the validation of non-extensive effects in
high-energy physics as well as in astrophysics bediaga ; lavagno2000 ;
lavagno2007 ; beck ; wilk2000 ; lavagno2001 ; biro2005 ; biro2009 ;
lavagno2002 ; alberico2008 ; quarati ; carvalho ; livadiotis ; leubner . More
details and diverse applications of Tsallis statistics can be found in
tsallis_2 .

The purpose of this paper is to study a strongly interacting system whose
dynamics are governed by non-extensive statistics. We investigate the
difference between thermodynamic properties of quark matter when Tsallis
statistics is applied and compare it with the BG statistics. In addition, we
have also studied the phase transition of quark matter within the realm of
non-extensive statistics. For this, we have generalized the PCQMF model to its
non-extensive version. We explore situations where both chemical potential and
temperature are non zero which stipulate the effect of Tsallis statistics on
the whole ($T-\mu_{B}$) plane in the phase diagram. Other models, such as the
QHD model, have also been generalized to their non-extensive version to study
the properties of nuclear matter periera2007 . Non-extensive statistical
effects on the properties of quark matter have also been studied within a
generalized NJL model rozy2009 . Furthermore, a generalized version of the MIT
bag model has also been used to study the QCD phase diagram and stellar matter
properties cardoso . In shen , a non-extensive linear sigma model is used to
study the chiral phase transition. Recently, the PNJL model has also been
generalized to its non-extensive version to study the thermodynamic properties
and transport coefficients zhao2020 , phase transition and critical exponents
of QCD matter zhao2021  as well as baryon number fluctuations zhao2023 .
Moreover, non-extensive statistics have also been used to study fluctuations
and correlations in high-energy nuclear collisions alberico2000 . In Ref.
drago , non-extensive effects have been considered to study the relativistic
nuclear equation of state. Ref. lavagno2010  included non-extensive effects to
study the hadron to QGP transition. The bulk properties of protoneutron stars
have also been investigated using the non-extensive formalism lavagno2011 .
Non-extensive thermodynamics of nuclear matter, as well as a black hole, has
been studied in megias . Moreover, in Ref. tiwari  non-extensive effects have
been used to study the transport properties of hadronic matter.

This paper is organized as follows: In Sec. II, we introduce the non-extensive
version of the Polyakov chiral SU (3) quark mean field model ($q$-PCQMF). The
impact of the non-extensive parameter, $q$, on various thermodynamic
properties of quark matter is discussed in Sec. III.1. In Sec. III.2, we
analyze in detail the influence of the parameter $q$ on the chiral and
deconfinement phase transition at finite quark chemical potential and finite
temperature. Finally, in Sec. IV, we provide a brief summary of our work.

##  II  Polyakov chiral SU(3) quark mean field model

Before presenting the $q$-PCQMF model, let us introduce the main concepts of
the PCQMF model. This model uses quarks and mesons as degrees of freedom. At
finite temperature and density, the interactions between quarks in this model
are described by the exchange of scalar meson fields $\sigma$, $\zeta$, and
$\delta$ and vector meson fields $\omega$, $\rho$ and $\phi$. The polyakov
fields $\Phi$ and $\bar{\Phi}$ are introduced in the model to study the
features of deconfinement phase transition. The attractive part of interaction
is represented by the scalar meson fields while the repulsive interaction is
through the vector meson fields. The non-strange scalar isoscalar meson field,
$\sigma$, corresponds to the scalar meson of mass 417 MeV wang2001  and has
($u\bar{u}/d\bar{d}$) quark as its content. The scalar meson interacts with
the strange quark through a strange scalar isoscalar field, $\zeta$ with
($s\bar{s}$) as its content and a mass of 1170 MeV wang2001 . The scalar
isovector field, $\delta$ is introduced in the model to study the effect of
isospin asymmetry and has ($u\bar{u}-d\bar{d}$) as its content. Furthermore,
the model introduces the broken scale invariance through a scalar dilaton
field, $\chi$ also known as glueball field schechter ; gomm ; heide ; carter ;
ko .

The effective Lagrangian of 2 + 1 flavour and 3 colour PCQMF model is as
follows

${\cal L}_{{\rm PCQMF}}\,=\,{\cal L}_{q0}\,+\,{\cal L}_{qm}\,+\,{\cal
L}_{M}\,+\,{\cal L}_{\Delta m}\,+\,{\cal
L}_{h}\,-U(\Phi(\vec{x}),\bar{\Phi}(\vec{x}),T).$ (3)

In the above equation ${\cal
L}_{q0}=\bar{\psi}\,i\gamma^{\mu}\partial_{\mu}\psi$, represents the free part
of massless quarks. The interactions between quark and meson is represented by

$\displaystyle{\cal
L}_{qm}=g_{s}\left(\bar{\psi}_{L}M\psi_{R}+\bar{\psi}_{R}M^{+}\psi_{L}\right)-g_{v}\left(\bar{\psi}_{L}\gamma^{\mu}l_{\mu}\psi_{L}+\bar{\psi}_{R}\gamma^{\mu}r_{\mu}\psi_{R}\right),$
(4)

where $\psi=\left(u,d,s\right)$ and $g_{v}$$(g_{s}$) are vector(scalar)
coupling constants. The spin-0 scalar ($\Sigma$) and pseudoscalar ($\Pi$)
mesons are expressed as papa1999

$M(M^{\dagger})=\Sigma\pm i\Pi=\frac{1}{\sqrt{2}}\sum_{a=0}^{8}(\sigma^{a}\pm
i\pi^{a})\lambda^{a},$ (5)

where $\sigma^{a}(\pi^{a})$ represent the scalar (pseudoscalar) meson nonets
and $\lambda^{a}$ are the Gell-Mann matrices with
$\lambda^{0}=\sqrt{\frac{2}{3}}\mathds{1}$. Similarly, spin-1 vector
($V_{\mu}$) and pseudovector ($A_{\mu}$) mesons are expressed as papa1999

$l_{\mu}(r_{\mu})=\frac{1}{2}(V_{\mu}\pm
A_{\mu})=\frac{1}{2\sqrt{2}}\sum_{a=0}^{8}(v^{a}_{\mu}\pm
a^{a}_{\mu})\lambda^{a},$ (6)

with $v^{a}_{\mu}$ ($a^{a}_{\mu}$) being vector (pseudovector) mesons nonets.

The meson interactions are described by the term ${\cal L}_{M}={\cal
L}_{\Sigma\Sigma}+{\cal L}_{VV}+{\cal L}_{SB}$. The self-interaction term for
the scalar meson in the mean field approximation can be expressed as

$\displaystyle{\cal L}_{\Sigma\Sigma}=$
$\displaystyle-\frac{1}{2}\,k_{0}\chi^{2}\left(\sigma^{2}+\zeta^{2}+\delta^{2}\right)+k_{1}\left(\sigma^{2}+\zeta^{2}+\delta^{2}\right)^{2}+k_{2}\left(\frac{\sigma^{4}}{2}+\frac{\delta^{4}}{2}+3\sigma^{2}\delta^{2}+\zeta^{4}\right)$
(7) $\displaystyle+k_{3}\chi\left(\sigma^{2}-\delta^{2}\right)\zeta-
k_{4}\chi^{4}-\frac{1}{4}\chi^{4}{\rm
ln}\frac{\chi^{4}}{\chi_{0}^{4}}+\frac{d}{3}\chi^{4}{\rm
ln}\left(\left(\frac{\left(\sigma^{2}-\delta^{2}\right)\zeta}{\sigma_{0}^{2}\zeta_{0}}\right)\left(\frac{\chi^{3}}{\chi_{0}^{3}}\right)\right),$

where $d=6/33$, $\sigma_{0}$ and $\zeta_{0}$ corresponds to the vacuum value
of $\sigma$ and $\zeta$ fields, respectively. These are expressed as
$\sigma_{0}=-f_{\pi}=93$ MeV and
$\zeta_{0}=\frac{1}{\sqrt{2}}(f_{\pi}-2f_{K})=115$ MeV, where $f_{\pi}$ is the
pion decay constant and $f_{K}$ is the kaon decay constant.

The self-interaction term for the vector meson is written as

${\cal
L}_{VV}=\frac{1}{2}\,\frac{\chi^{2}}{\chi_{0}^{2}}\left(m_{\omega}^{2}\omega^{2}+m_{\rho}^{2}\rho^{2}+m_{\phi}^{2}\phi^{2}\right)+g_{4}\left(\omega^{4}+6\omega^{2}\rho^{2}+\rho^{4}+2\phi^{4}\right).$
(8)

where the density-dependence of the mass of vector meson can be written as
wang2003

$\displaystyle
m^{2}_{\omega}=m_{\rho}^{2}=\frac{m_{v}^{2}}{1-\frac{1}{2}\mu\sigma^{2}}\,,\hskip
14.22636pt{\rm and}\hskip
14.22636ptm^{2}_{\phi}=\frac{m_{v}^{2}}{1-\mu\zeta^{2}}\,.$ (9)

The parameters $m_{v}$ = 673.6 MeV and density parameter $\mu$ = 2.34 fm2 are
taken to yield ${m_{\phi}}$ = 1020 MeV and ${m_{\omega}}$ = 783 MeV.

The masses of pseudo-scalar mesons arises due to the term wang2002 ; wang2004
; papa1999

${\cal
L}_{SB}=-\frac{\chi^{2}}{\chi_{0}^{2}}\left[m_{\pi}^{2}f_{\pi}\sigma+\left(\sqrt{2}m_{K}^{2}f_{K}-\frac{m_{\pi}^{2}}{\sqrt{2}}f_{\pi}\right)\zeta\right].$
(10)

The exact mass of $s$ quark is generated through ${\cal L}_{\Delta m}=-\Delta
m_{s}\bar{\psi}S\psi$, where $\Delta m_{s}=29$ MeV and
$S\,=\,\frac{1}{3}\,\left(I-\lambda_{8}\sqrt{3}\right)={\rm diag}(0,0,1)$.

In order to study the dependence of $T$ and $\mu_{B}$ in this model, we
investigate the grand canonical potential density which can be expressed in
the mean field approximation as

$\Omega_{\rm{PCQMF}}=\mathcal{U}(\Phi,\bar{\Phi},T)-{\cal L}_{M}-{\cal
V}_{vac}+\sum_{i=u,d,s}\frac{-\gamma_{i}k_{B}T}{(2\pi)^{3}}\int_{0}^{\infty}d^{3}k\left\\{{\rm
ln}F^{-}+{\rm ln}F^{+}\right\\},$ (11)

where, $\gamma_{i}$ = 2 is the spin degeneracy factor, and

$\displaystyle F^{-}=$ $\displaystyle 1+e^{-3E^{-}}+3\Phi
e^{-E^{-}}+3\bar{\Phi}e^{-2E^{-}},$ (12) $\displaystyle F^{+}=$ $\displaystyle
1+e^{-3E^{+}}+3\bar{\Phi}e^{-E^{+}}+3\Phi e^{-2E^{+}}.$ (13)

In the above equation, $E^{-}=(E_{i}^{*}(k)-{\nu_{i}}^{*})/k_{B}T$ and
$E^{+}=(E_{i}^{*}(k)+{\nu_{i}}^{*})/k_{B}T$ in which the effective single
particle energy of quarks is given by $E_{i}^{*}(k)=\sqrt{m_{i}^{*2}+k^{2}}$.
In Eq. (11) the term, ${\cal V}_{vac}$ is subtracted for the sake of getting
zero vacuum energy. The in-medium quark chemical potential ${\nu_{i}}^{*}$ is
linked to the chemical potential $\nu_{i}$ in free space by

${\nu_{i}}^{*}=\nu_{i}-g_{\omega}^{i}\omega-g_{\phi}^{i}\phi-
g_{\rho}^{i}\rho,$ (14)

where, $g^{i}_{\omega}$, $g^{i}_{\phi}$ and $g^{i}_{\rho}$ are the coupling
constants of vector meson fields with different quarks. The effective mass of
quark, ${m_{i}}^{*}$, is given by

${m_{i}}^{*}=-g_{\sigma}^{i}\sigma-g_{\zeta}^{i}\zeta-
g_{\delta}^{i}\delta+{m_{i0}},$ (15)

where $g_{\sigma}^{i}$, $g_{\zeta}^{i}$ and $g_{\delta}^{i}$ describe the
coupling strength of scalar fields with different quarks. The values of
$g_{\sigma}^{i}$, $g_{\zeta}^{i}$ and ${m_{i0}}$ are determined by fitting the
vacuum masses of constituent quarks, $m_{u},m_{d},m_{s}$ = 313 MeV, 313 MeV,
490 MeV, respectively wang2003 .

For the Polyakov loop effective potential, we have considered the commonly
used logarithmic form given by fukushima ; costa ; roessner ,

$\displaystyle\frac{U(\Phi,\bar{\Phi},T)}{T^{4}}$ $\displaystyle=$
$\displaystyle-\frac{a(T)}{2}\bar{\Phi}\Phi+b(T)\mathrm{ln}\big{[}1-6\bar{\Phi}\Phi+4(\bar{\Phi}^{3}+\Phi^{3})-3(\bar{\Phi}\Phi)^{2}\big{]},$
(16)

with the temperature dependent coefficients costa ; roessner2008 :

$a(T)=a_{0}+a_{1}\bigg{(}\frac{T_{0}}{T}\bigg{)}+a_{2}\bigg{(}\frac{T_{0}}{T}\bigg{)}^{2},\
\ b(T)=b_{3}\bigg{(}\frac{T_{0}}{T}\bigg{)}^{3}.$ (17)

The corresponding parameters $a_{0}$, $a_{1}$, $a_{2}$ and $b_{3}$ are given
in Table 1.

$a_{0}$ | $a_{1}$ | $a_{2}$ | $b_{3}$
---|---|---|---
1.81 | -2.47 | 15.2 | -1.75
Table 1: Parameters in Polyakov effective potential

In the pure gauge sector at vanishing chemical potential, the parameter
$T_{0}=270$ MeV fukugita . $T_{0}$ is rescaled from 270 MeV to 200 MeV when
fermionic fields are included to compare the model results with existing
lattice results manisha . Finally, the solutions of the scalar fields
$\sigma$, $\zeta$ and $\delta$, the dilaton field $\chi$, the vector fields
$\omega$, $\rho$ and $\phi$ and the Polyakov field $\Phi$ and its conjugate
$\bar{\Phi}$ are obtained by minimizing $\Omega$ with respect to these fields,
$i.e.$,

$\frac{\partial\Omega}{\partial\sigma}=\frac{\partial\Omega}{\partial\zeta}=\frac{\partial\Omega}{\partial\delta}=\frac{\partial\Omega}{\partial\chi}=\frac{\partial\Omega}{\partial\omega}=\frac{\partial\Omega}{\partial\rho}=\frac{\partial\Omega}{\partial\phi}=\frac{\partial\Omega}{\partial\Phi}=\frac{\partial\Omega}{\partial\bar{\Phi}}=0.$
(18)

The model parameters which are summarized in Table 2. can be estimated using
the masses of $\pi,K$ mesons, average masses of $\eta$ and
$\eta^{{}^{\prime}}$ as well as the vacuum masses of $\sigma$, $\zeta$ and
$\chi$ mesons.

$k_{0}$ | $k_{1}$ | $k_{2}$ | $k_{3}$ | $k_{4}$ | $g_{s}$ | $\rm{g_{v}}$ | $\rm{g_{4}}$ | $h_{1}$ | $h_{2}$
---|---|---|---|---|---|---|---|---|---
4.94 | 2.12 | -10.16 | -5.38 | -0.06 | 4.76 | 10.92 | 37.5 | 0 | 0
$\sigma_{0}$ (MeV) |  $\zeta_{0}$(MeV) |  $\chi_{0}$(MeV) |  $m_{\pi}$(MeV) |  $f_{\pi}$(MeV) |  $m_{K}$(MeV) |  $f_{K}$(MeV) |  $m_{\omega}$(MeV) |  $m_{\phi}$(MeV) |  $m_{\rho}$( MeV)
-93 | -96.87 | 254.6 | 139 | 93 | 496 | 115 | 783 | 1020 | 783
$g_{\sigma}^{u}$ | $g_{\sigma}^{d}$ | $g_{\sigma}^{s}$ | $g_{\zeta}^{u}$ | $g_{\zeta}^{d}$ | $g_{\zeta}^{s}$ | $g_{\delta}^{u}$ | $g_{\delta}^{d}$ | $g_{\delta}^{s}$ |  $\rho_{0}$(fm-3)
3.36 | 3.36 | 0 | 0 | 0 | 4.76 | 3.36 | -3.36 | 0 | 0.15
$g^{u}_{\omega}$ | $g^{d}_{\omega}$ | $g^{s}_{\omega}$ | $g^{u}_{\phi}$ | $g^{d}_{\phi}$ | $g^{s}_{\phi}$ | $g^{u}_{\rho}$ | $g^{d}_{\rho}$ | $g^{s}_{\rho}$ | $d$
3.86 | 3.86 | 0 | 0 | 0 | 5.46 | 3.86 | -3.86 | 0 | 0.18
Table 2: The list of parameters used in the present work.

###  II.1  q-PCQMF model

The PCQMF model assumes the additivity of some thermodynamical quantities,
like, entropy. For systems where long-range correlations and fluctuations are
important, e.g., under phase transition, this is a very strong approximation.
One way to take the nonadditivity of interacting systems into account is to
consider a micro-canonical ensemble of Hamiltonian systems gross2004 ;
gross2006 . An alternative approach is to make use of Tsallis statistics
tsallis . As discussed in Sec. (I), resigning from the assumption of
additivity and describing theoretical models based on the non-extensive
approach appears to be a promising possibility. In this spirit, we have
extended the PCQMF model to its non-extensive version by using Tsallis
statistics instead of BG statistics, which implies the replacements described
in Eq.(1) and (2). However, we have taken the following simplifications:
Firstly, Polyakov loop potential experience non-extensive effects only
implicitly and hence remain unchanged. Secondly, similar to the case of
finite-size effects where the volume $V$ is treated on the same footing as $T$
and $\mu$ nisha , the usual parameters of the PCQMF model are left undisturbed
and the $q$ parameter is treated on the same footing as $T$ and $\mu$.

Hence, within the $q$-PCQMF model, the thermodynamic potential density becomes

$\Omega_{q}=\mathcal{U}(\Phi,\bar{\Phi},T)-{\cal L}_{M}-{\cal
V}_{vac}+\sum_{i=u,d,s}\frac{-\gamma_{i}k_{B}T}{(2\pi)^{3}}\int_{0}^{\infty}d^{3}k\left\\{{\rm
ln_{q}}F_{q}^{-}+{\rm ln_{q}}F^{+}_{q}\right\\},$ (19)

where

$\displaystyle F_{q}^{-}=$ $\displaystyle 1+e_{q}(-3E^{-})+3\Phi
e_{q}(-E^{-})+3\bar{\Phi}e_{q}(-2E^{-}),$ (20) $\displaystyle F_{q}^{+}=$
$\displaystyle 1+e_{q}(-3E^{+})+3\bar{\Phi}e_{q}(-E^{+})+3\Phi
e_{q}(-2E^{+}).$ (21)

The coupled equations obtained by minimizing $\Omega_{q}$ in Eq. (19) are
given in the appendix. The $q$ versions of vector density, $\rho_{q,i}$, and
scalar density, $\rho_{q,i}^{s}$, of quarks are written as

$\rho_{q,i}=\gamma_{i}N_{c}\int\frac{d^{3}k}{(2\pi)^{3}}\Big{(}f_{q,i}(k)-\bar{f}_{q,i}(k)\Big{)},$
(22)

and

$\rho_{q,i}^{s}=\gamma_{i}N_{c}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{m_{i}^{*}}{E^{\ast}_{i}(k)}\Big{(}f_{q,i}(k)+\bar{f}_{q,i}(k)\Big{)},$
(23)

respectively, where $f_{q,i}(k)$ and $\bar{f}_{q,i}(k)$ are the $q-$modified
Fermi-distribution functions for quark and anti-quark at finite temperature
and are given by,

$f_{q,i}(k)=\frac{\Phi
e^{q}_{q}(-E^{-})+2\bar{\Phi}e^{q}_{q}(-2E^{-})+e^{q}_{q}(-3E^{-})}{[1+3\Phi
e_{q}(-E^{-})+3\bar{\Phi}e_{q}(-2E^{-})+e_{q}(-3E^{-})]^{q}},$ (24)

$\bar{f}_{q,i}(k)=\frac{\bar{\Phi}e^{q}_{q}(-E^{+})+2\Phi
e^{q}_{q}(-2E^{+})+e^{q}_{q}(-3E^{+})}{[1+3\bar{\Phi}e_{q}(-E^{+})+3\Phi
e_{q}(-2E^{+})+e_{q}(-3E^{+})]^{q}}.$ (25)

It is important to note that in the limit $q\rightarrow 1,$ the standard
Fermi-distribution functions are recovered and hence we return back to the
standard (extensive) PCQMF model.

The following restriction must hold for $e_{q}(x)$ to always be a non-negative
real function:

$[1+(1-q)x]\geq 0$ (26)

In this paper, we consider only $q>1$. This is due to the fact that it is
discovered that the typical value of the non-extensivity parameter $q$ for
high-energy collisions is $1\leq q\leq 1.2$ marques ; li2013 ; cleymans2013 ;
azmi . Furthermore, for $q>1$, $q-1$ describes intrinsic fluctuations of
temperature in the system wilk2000 ; biro2005 , whereas, the interpretation of
$q-1$ for $q<1$ is inconsistent kodama ; garcia . To comply with Eq.(26), for
$q>1$, we can use the following Tsallis cutoff prescription tewel

$e_{q}(x)=\begin{cases}[1+(1-q)x]^{\frac{1}{1-q}},\quad\text{for}\quad x\leq
0,\\\ [1+(q-1)x]^{\frac{1}{q-1}},\quad\text{for \quad}x>0\quad.\end{cases}$
(27)

Also, let us observe that as $T\rightarrow 0$, as long as $q>1$,
$\Omega_{q}\rightarrow\Omega$ and the non-extensive effects vanish. This means
that in studies where in the nuclear scale $T\approx 0$ (in the interior of
neutron stars) we expect no non-extensive signature. Whereas, when $T$ is in
MeVs, (heavy-ion collision experiments), one expects non-extensive effects to
play an important role in the thermodynamic quantities. To study such effects,
we apply the modified thermodynamic potential density Eq.(19) in Eq.(18) to
calculate the in-medium $q$-value of the vector and scalar fields. The isospin
asymmetry can be incorporated through the definition
$\eta=\frac{(\rho_{d}-\rho_{u})}{(\rho_{d}+\rho_{u})/3}$. The total baryon
density is given by $\rho_{B}=\frac{1}{3}(\rho_{u}+\rho_{d}+\rho_{s})$. The
various chemical potentials are defined through
$\mu_{B}=\frac{3}{2}(\mu_{u}+\mu_{d})$, $\mu_{I}=\frac{1}{2}(\mu_{u}-\mu_{d})$
and $\mu_{S}=\frac{1}{2}(\mu_{u}+\mu_{d}-2\mu_{s})$, correspond to baryon,
isospin and strangeness chemical potential, respectively.

##  III  Results and Discussions

In this section, we present various results on the thermodynamical properties
of quark matter within $q-$extended Polyakov chiral SU(3) quark mean-field
model ($q$-PCQMF) presented in Sec. (II.1). In the $q$-PCQMF model, the effect
of the $q$ parameter comes into the picture through the thermodynamic
potential density, $\Omega_{q}$, which depends upon the $q-$ dependent scalar
density, $\rho_{i}^{s}$ and vector density, $\rho_{i}$ of constituent quarks,
which in turn depend upon the scalar fields ($\sigma$, $\zeta$, $\delta$,
$\chi$), the vector fields ($\omega$, $\rho$, $\phi$), and the Polyakov loop
fields, ($\Phi$, $\bar{\Phi}$). As said before, these fields are calculated by
solving the coupled system of non-linear equations. In Sect. III.1, the in-
medium properties of the scalar and the Polyakov loop fields will be
discussed. These fields will be used as input to understand the effect of $q$
on different thermodynamic properties of quark matter. We will discuss the
behaviour of the susceptibility of these fields and study the $(T-\mu)$ phase
diagram and the influence of the $q$ parameter on these quantities in Sect.
III.2.

###  III.1  In-medium fields and thermodynamic properties

In Fig. (1), the variations of the scalar fields $\sigma$ and $\zeta$ as well
as the Polyakov fields $\Phi$ and $\bar{\Phi}$ are shown as a function of
temperature, $T$, for $q$ = 1, 1.05, and 1.10 at baryon chemical potential,
$\mu_{B}=400$ MeV, isospin chemical potential, $\mu_{I}=40$ MeV and
strangeness chemical potential, $\mu_{S}=125$ MeV. We observe that the
magnitude of $\sigma$ and $\zeta$ decrease with an increase in temperature.
The temperature at which the magnitude of scalar fields begins to decrease
rapidly is referred to as pseudo-critical temperature, $T_{p}$. This decrease
in the magnitude of scalar fields may represent the restoration of chiral
symmetry at high temperatures. Here, we can see that as $q$ increases, the
value of $T_{p}$ decreases. As the temperature is increased the Fermi
distribution function in Eqs. (24) and (25) decreases resulting in a decrease
in the magnitude of the scalar fields. At higher values of $q$, there is a
considerable decrease in the magnitude of the fields even at lower
temperatures. This may signify that the transition temperature is shifted
towards a lower temperature value.

In the mean-field approximation, the order parameters to study the
deconfinement are the Polyakov loop fields, $\Phi$ and $\bar{\Phi}$. Variation
in the magnitude of $\Phi$ and $\bar{\Phi}$ gives information about the
deconfinement phase transition. As can be seen in Fig. 1 (c) and (d), at lower
temperatures, the value of both $\Phi$ and $\bar{\Phi}$ is approximately zero,
indicating confined quarks within hadrons. With an increase in temperature,
$\Phi$ and $\bar{\Phi}$ increase and reaches a critical temperature, where it
jumps to a non-zero value, indicating the deconfinement phase transition.
Similar to the scalar fields, we find that as non-extensivity increases
(higher values of $q$), this increase in $\Phi$ and $\bar{\Phi}$ occurs
earlier. This may indicate a decrease in the deconfinement temperature for
$q>1$.

Figure 1: Scalar fields $\sigma$, $\zeta$ and Polyakov fields $\Phi$,
$\bar{\Phi}$ plotted as a function of temperature T, for $q$ = 1, 1.05, and
1.10 at $\mu_{B}=400$ MeV, $\mu_{I}=40$ MeV, $\mu_{S}=125$ MeV.

One can calculate various thermodynamic properties such as energy density,
pressure density, entropy density, and trace anomaly of quark matter with the
help of the thermodynamics potential density, $\Omega_{q}$ in Eq. 19. The
pressure is given by

$p_{q}(T)=-\Omega_{q}(T).$ (28)

The entropy density and energy density are given by

$s_{q}(T)=-\frac{\partial\Omega_{q}(T)}{\partial
T},~{}~{}~{}~{}~{}~{}\epsilon_{q}(T)=\Omega_{q}(T)+\sum_{i=u,d,s}{\nu_{i}}^{*}\rho_{i}+Ts_{q}(T),$
(29)

respectively.

Fig. 2 shows the behaviour of $\epsilon_{q}/T^{4}$, $p_{q}/T^{4}$,
$s_{q}/T^{3}$ and $(\epsilon_{q}-3p_{q})/T^{4}$ at $\mu_{B}=400$ MeV,
$\mu_{I}=40$ MeV, $\mu_{S}=125$ MeV for $q=1,1.05$, and $1.10$. All these
quantities show smooth behaviour implying a crossover phase transition. We see
a sharp increase in the value of these quantities near the transition
temperature. It is interesting to note that when $q=1$, these thermodynamic
quantities all tend to their Stefan Boltzmann (SB) limit at higher
temperatures. However, with an increase in $q$, they increase rapidly and
exceed their corresponding SB limit. For $\epsilon_{q}/T^{4}$ at $T=300$ MeV,
the value at $q=1$ is 14.39 (very close to the SB limit, 15.627) but as $q$ is
increased to 1.1, its value increases by $75\%$ to 25.24. For $s_{q}/T^{3}$
and $p_{q}/T^{4}$, these values are $80\%$ and $105\%$, respectively. This is
the anticipated result in the $q-PCQMF$ model, because of the usage of Tsallis
statistics instead of BG statistics. In Tsallis statistics, $q-1$ is
interpreted as a deviation from the BG statistics biro20052 ; wilk2000 . The
higher the value of $q-1$, the more the deviation from the BG statistics and
correspondingly, the higher the value of the SB limit for quark matter. This
is because when $T\rightarrow\infty$, the modified thermodynamic potential
density doesn’t return to its normal value, i.e., $\Omega_{q}\neq\Omega$. This
means that for a system whose dynamics is described by Tsallis statistics, the
high-temperature limit of various thermodynamic quantity is no longer the SB
limit but a $q$-dependent Tsallis limit. This is similar to the results
obtained in the $q$-NJL model zhao2020 .

Some other relevant quantities in relativistic heavy-ion collision are the
speed of sound squared at constant entropy $c_{sq}^{2}$ and specific heat at
constant volume $c_{vq}$, which are defined as

$c_{sq}^{2}=\left(\frac{\partial
p_{q}}{\partial\epsilon_{q}}\right)_{s_{q}}=\frac{s_{q}}{c_{vq}},\hskip
30.0ptc_{vq}=\left(\frac{\partial\epsilon_{q}}{T}\right)_{V}.$ (30)

From Fig. 3 (a), we can see that the value of $c_{v}/T^{3}$ shows a similar
trend to the other thermodynamic quantities discussed above. At $q=1$, its
value increases with an increase in temperature reaches a maximum near the
pseudo-critical temperature and then tends to the usual SB limit, 63.14.
Besides, for higher values of $q$, it also reaches a new $q$-dependent Tsallis
limit. In Fig. 3 (b), at $q=1$, the value of $c_{s}^{2}$ dips near the pseudo-
critical temperature but then approaches the SB value, 0.33. However, as $q$
increases, the dip disappears and the criticality vanishes. It is interesting
to note that at high temperatures, both $s_{q}/T^{3}$ and $c_{vq}/T^{3}$ show
similar growth. Therefore, non-extensive effects such as high-temperature
limit are not observed in
$c_{sq}^{2}=\frac{s_{q}}{T^{3}}/\frac{c_{vq}}{T^{3}}$. Similar results are
obtained in Ref. khuntia .

Figure 2: Scaled energy density, $\epsilon_{q}/T^{4}$, pressure density,
$p_{q}/T^{4}$, entropy density, $s_{q}/T^{3}$ and trace anomaly,
$(\epsilon_{q}-3p_{q})/T^{4}$ as a function of temperature $T$, for $q$ = 1,
1.05, and 1.10 at $\mu_{B}=400$ MeV, $\mu_{I}=40$ MeV, $\mu_{S}=125$ MeV.
Figure 3: Scaled specific heat, $c_{vq}/T^{3}$ and the square of the speed of
sound, $c_{sq}^{2}$ as a function of temperature $T$, for $q$ = 1, 1.05, and
1.10 at $\mu_{B}=400$ MeV, $\mu_{I}=40$ MeV, $\mu_{S}=125$ MeV.

###  III.2  $T-\mu$ phase diagram

In this subsection, we study the $T-\mu$ phase diagram within the $q$-PCQMF
model. We analyse three phase transitions in the present work: the chiral
phase transition from the chirally broken phase to the chiral restoration
phase for up and down quarks, the chiral phase transition for the strange
quark, and the deconfinement phase transition. The effect of the $q$ parameter
on the nature of these phase transitions is studied with the help of the
susceptibilities, defined as

$\chi_{\sigma_{x}}=\left(\frac{\partial\sigma_{x}}{\partial
T}\right)_{\mu},\hskip
25.0pt\chi_{\sigma_{y}}=\left(\frac{\partial\sigma_{y}}{\partial
T}\right)_{\mu},\hskip 25.0pt\chi_{\Phi}=\left(\frac{\partial\Phi}{\partial
T}\right)_{\mu},\hskip
25.0pt\chi_{\bar{\Phi}}=\left(\frac{\partial\bar{\Phi}}{\partial
T}\right)_{\mu},$ (31)

where $\sigma_{x}$, for $u$ and $d$ quarks, and $\sigma_{y}$, for $s$ quark,
are the scalar quark condensates which can be expressed in terms of scalar
fields $\sigma,\zeta,\delta,$ and $\chi$ as manisha ,

$\sigma_{x}=\sigma_{u,d}=\frac{1}{m_{u,d}}\left(\frac{\chi}{\chi_{0}}\right)^{2}\Bigg{[}\frac{1}{2}m^{2}_{\pi}f_{\pi}(\sigma\pm\delta)\Bigg{]},$
(32)
$\sigma_{y}=\frac{1}{m_{s}}\left(\frac{\chi}{\chi_{0}}\right)^{2}\Bigg{[}\left(\sqrt{2}m^{2}_{k}f_{k}-\frac{1}{\sqrt{2}}m^{2}_{\pi}f_{\pi}\right)\zeta\Bigg{]}.$
(33)

Figure 4: Susceptibilities of $\sigma_{x}$, $\sigma_{y}$, $\Phi$, and
$\bar{\Phi}$ at $q$ = 1, 1.05, and 1.10 for $\mu_{B}=400$ MeV, $\mu_{I}=40$
MeV, $\mu_{S}=125$ MeV.

From Fig. (1) it is evident that there is a smooth transition of the order
parameters (scalar fields ($\sigma$ and $\zeta$) and Polyakov loop ($\Phi$ and
$\bar{\Phi}$)). We do however see a change in the behaviour of these order
parameters as we increase the temperature. This is evidence of a crossover
phase transition and we can study the susceptibilities defined in Eq. (31) to
figure out the position of the pseudo-critical point $T_{p}$. The peak value
of the susceptibilities $\chi_{\sigma_{x}},\chi_{\sigma_{y}}$ corresponds to
the pseudo-critical temperatures $T^{q}_{\chi}$, $T_{\chi}^{s}$ for chiral
phase transition of up/down and strange quarks, respectively. The pseudo-
critical temperature for the deconfinement phase transition $T_{d}$ is
obtained through $\chi_{\Phi}$. Fig. (4) shows the variation of
susceptibilities $\chi_{\sigma_{x}}$, for up and down quarks,
$\chi_{\sigma_{y}}$, for the strange quark, and the Polyakov loop,
$\chi_{\Phi}$ and $\chi_{\bar{\Phi}}$ at $\mu_{B}=400,\mu_{I}=40,\mu_{S}=125$
MeV for $q=1,1.05$, and $1.10$. For $q$ = 1, $\chi_{\sigma_{x}}$ show a peak
at $T\approx 160$ MeV, whereas, $\chi_{\sigma_{y}}$ has two peaks, one at
$T\approx 160$ MeV and another around $T\approx 247$ MeV. Furthermore, only a
single peak which coincides with the peak of $\chi_{\sigma_{x}}$ at $T\approx
160$ MeV is observed for $\chi_{\Phi}$ and $\chi_{\bar{\Phi}}$ at $q=1$. The
peak of the susceptibility gives us the pseudo-critical temperature. In case
of two peaks in susceptibilities, the criteria for deciding the pseudo-
critical temperature for chiral phase transition is given by
$\sigma_{x,y}(T)/\sigma_{x,y}(T=0)<1/2$, and for the deconfinement phase
transition is $\Phi(T)/\Phi(T\rightarrow\infty)>1/2$ mao . It is evident from
Fig. (4) that the value of these pseudo-critical temperatures decreases with
increasing $q$ for both the chiral phase transition $T_{\chi}^{q}$ and the
deconfinement phase transition $T_{d}$. It is also interesting to note that
the chiral phase transition for the strange quark occurs at a much higher
temperature as evidenced by the second peak at $T_{\chi}^{s}\approx 247$ MeV,
whereas for lighter up and down quark it occurs at $T_{\chi}^{q}\approx 160$
MeV.

Figure 5: Quark number density $\rho_{u}$ and susceptibility $\chi_{u}$ of up
quark as a function of quark chemical potential at T = 70, 80, 90, and 100 MeV
for $q=1$. Figure 6: Quark number density $\rho_{u}$ and susceptibility
$\chi_{u}$ of up quark as a function of quark chemical potential at T = 70,
80, 90, and 100 MeV for $q=1.05$. Figure 7: Quark number density $\rho_{u}$
and susceptibility $\chi_{u}$ of up quark as a function of quark chemical
potential at T = 70, 80, 90, and 100 MeV for $q=1.10$.

The nature of the phase transition at a higher value of quark chemical
potential is explored using quark number susceptibility. The chiral phase
transition changes from a crossover to a first-order phase transition at
higher values of $\mu_{q}$ and this can be seen by studying the discontinuity
in the quark number density, $\rho_{u}$, or divergence in the corresponding
susceptibility, $\chi_{u}$ defined as

$\chi_{u}=\left(\frac{\partial\rho_{u}}{\partial\mu}\right)_{T}.$ (34)

In Fig. 5 (a), quark number density is plotted as a function of quark chemical
potential for different values of temperature and zero vector interaction, at
$q$ = 1. It can be seen that when $T\leq T_{CEP}$, $\rho_{u}$ shows a
discontinuity, which can be seen in the corresponding susceptibility,
$\chi_{u}$, plotted in Fig. 5 (b), indicating the presence of first-order
phase transition. When the temperature rises above $T_{CEP}$, the
discontinuity in $\rho_{u}$ disappears suggesting the order of the phase
transition to be a crossover. This change in the nature of $\rho_{u}$ can be
used to determine the location of the CEP. For $q$ = 1, we observed the
position of the CEP at around $(T,\mu)$=(90 MeV, 286 MeV) manisha . When the
value of $q$ is increased to 1.05, the peak in the $\chi_{u}$ is shifted
towards higher $\mu$ and lower $T$ as shown in Fig. 6. For $q$ = 1.05, the
position of the CEP is around $(T,\mu)$=(84 MeV, 290 MeV). A similar trend is
shown in Fig. (7) for $q=1.10$, where we found the CEP at around $(T,\mu)$=(78
MeV, 293 MeV).

Fig. (8) shows the deconfinement phase transition at $\mu_{S}=\mu_{I}=0$ for
different $q$ values. This deconfinement phase transition remains a crossover
over the entire range of the quark chemical potential. As the value of $q$
increases, the phase transition temperature decreases to lower values. This
means that with increasing non-extensivity, quarks can become deconfined even
at lower temperatures. For $q$ = 1, at zero chemical potential, the value of
this temperature, $T_{d}$ is 164 MeV. This value decreases to $T_{d}$ = 153
MeV for $q$ = 1.05 and $T_{d}$ = 140 MeV for $q$ = 1.10. However, the decrease
is less significant at a higher value of $\mu$. For example, at $\mu=350$ MeV,
the values of $T_{d}$ are observed to be 116, 111, and 105 MeV. This signifies
that the impact of non-extensivity vanishes as we move towards lower
temperatures. In Fig. 9 (a), the chiral phase transition of up and down quark
is shown for the $q$-PCQMF model in the $T-\mu$ plane for different values of
$q$, at $\mu_{S}=\mu_{I}=0$. For plotting the phase boundary in the crossover
region, we have used the maximum of $\chi_{\sigma}$ for a given value of
$\mu$, whereas, for the first-order transition we have used the maximum of
$\chi_{u}$ for a given value of $T$. With an increase in non-extensivity of
the system, the regime of first-order phase transition becomes smaller and the
CEP is achieved at a higher value of $\mu$ and lower value of $T$. It is also
clear from Fig. 9 (a) that the crossover transition, which occurs at higher
temperatures is distinctly different but the first-order phase transition,
which occurs at lower temperatures is similar when $q>1$. This is again
evidence of the fact that non-extensivity vanishes at lower temperatures,
i.e., $\Omega_{q}\rightarrow\Omega$ when $T\rightarrow 0$. The chiral phase
transition of the s quark is shown in Fig. 9 (b), for $q$ = 1, 1.05, and 1.10.
Contrary to the chiral phase transition of the u and d quark, this phase
transition remains a crossover over the entire range of $(T-\mu)$ phase
diagram. Besides, this phase transition also occurs at lower temperature
values as $q$ becomes greater than one. However, this difference in
temperature at different $q$ values disappears towards lower temperatures or
higher chemical potentials indicating the vanishing of non-extensive effects.

Figure 8: Deconfinement phase transition of quark matter at $\mu_{S}$ =
$\mu_{I}$ = 0 MeV for $q$ = 1, 1,05, and 1.10. Figure 9: Chiral phase
transition: (a) for $u$ and $d$ quarks, (b) for $s$ quark, at $\mu_{S}$ =
$\mu_{I}$ = 0 MeV for $q$ = 1, 1,05, and 1.10.

##  IV  Summary and Conclusions

To summarize, we have extended the PCQMF model to its non-extensive version
and studied the properties of quark matter at finite temperature and density.
We have studied the thermodynamic properties and phase transitions with
Tsallis statistics and compared it to the usual BG statistics. The non-
extensivity is introduced into our model through the thermodynamic potential
density which modifies the form of the scalar and the vector densities of the
quark resulting in the modification of the scalar, vector as well as polyakov
loop fields. These modified fields are then further used to calculate various
thermodynamic quantities. We have found that the quantities such as
$\frac{\epsilon_{q}}{T^{4}},\frac{p_{q}}{T^{4}},\frac{s_{q}}{T^{3}},$ and
$\frac{(\epsilon_{q}-3p_{q})}{T^{4}}$ all tend towards a new $q$-related
Tsallis limit rather than their usual SB limit at high temperature for $q>1$.
We found that as $q$ increases, the criticality of $c_{v}/T^{3}$ and
$c_{s}^{2}$ gradually decreases. Also, the high-temperature limit of
$c_{s}^{2}$ is unaffected by the $q$ parameter due to a surprising
cancellation zhao2020 . Regarding the influence of the $q$ parameter on the
deconfinement phase transition, we found that this phase transition occurs at
a lower value of $T$ as the $q$ parameter is increased. However, the nature of
this phase transition is still a crossover, independent of $q$. Further, we
have also studied the influence of the $q$ parameter on the chiral phase
transition for the $u$ and $d$ quarks as well as for the s quark. We found
that for $u$ and $d$ quarks, the nature of this phase transition changes order
from a crossover at high temperatures to a first-order phase transition at
lower temperatures. Also, as $q$ increases the CEP is shifted towards a lower
$T$ and a higher $\mu$ value. However, for the s quark, the phase transition
remains a crossover but is shifted towards lower $T$ values as $q$ becomes
greater than one. In addition, we have found that the non-extensivity of the
system disappears as we approach $T\rightarrow 0$. In future, we aim to study
fluctuations and correlations of conserved charges on the nuclear and
hyperonic matter in the realm of non-extensivity.

## Appendix

The coupled equations are obtained by minimizing the thermodynamic potential
density, $\Omega_{q}$ with respect to the various fields of the $q$-PCQMF
model and are given as

$\displaystyle\frac{\partial\Omega_{q}}{\partial\sigma}=k_{0}\chi^{2}\sigma-4k_{1}\left(\sigma^{2}+\zeta^{2}+\delta^{2}\right)\sigma-2k_{2}\left(\sigma^{3}+3\sigma\delta^{2}\right)-2k_{3}\chi\sigma\zeta$
(35) $\displaystyle-$
$\displaystyle\frac{d}{3}\chi^{4}\bigg{(}\frac{2\sigma}{\sigma^{2}-\delta^{2}}\bigg{)}+\left(\frac{\chi}{\chi_{0}}\right)^{2}m_{\pi}^{2}f_{\pi}-\left(\frac{\chi}{\chi_{0}}\right)^{2}m_{\omega}\omega^{2}\frac{\partial
m_{\omega}}{\partial\sigma}$ $\displaystyle-$
$\displaystyle\left(\frac{\chi}{\chi_{0}}\right)^{2}m_{\rho}\rho^{2}\frac{\partial
m_{\rho}}{\partial\sigma}-\sum_{i=u,d}g_{\sigma}^{i}\rho_{q,i}^{s}=0,$
$\displaystyle\frac{\partial\Omega_{q}}{\partial\zeta}=k_{0}\chi^{2}\zeta-4k_{1}\left(\sigma^{2}+\zeta^{2}+\delta^{2}\right)\zeta-4k_{2}\zeta^{3}-k_{3}\chi\left(\sigma^{2}-\delta^{2}\right)-\frac{d}{3}\frac{\chi^{4}}{{\zeta}}$
(36) $\displaystyle+$
$\displaystyle\left(\frac{\chi}{\chi_{0}}\right)^{2}\left[\sqrt{2}m_{K}^{2}f_{K}-\frac{1}{\sqrt{2}}m_{\pi}^{2}f_{\pi}\right]-\left(\frac{\chi}{\chi_{0}}\right)^{2}m_{\phi}\phi^{2}\frac{\partial
m_{\phi}}{\partial\zeta}-\sum_{i=s}g_{\zeta}^{i}\rho_{q,i}^{s}$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\frac{\partial\Omega_{q}}{\partial\delta}=k_{0}\chi^{2}\delta-4k_{1}\left(\sigma^{2}+\zeta^{2}+\delta^{2}\right)\delta-2k_{2}\left(\delta^{3}+3\sigma^{2}\delta\right)+\mathrm{2k_{3}\chi\delta\zeta}$
(37) $\displaystyle+$
$\displaystyle\frac{2}{3}d\chi^{4}\left(\frac{\delta}{\sigma^{2}-\delta^{2}}\right)-\sum_{i=u,d}g_{\delta}^{i}\rho_{q,i}^{s}=0,$
$\displaystyle\frac{\partial\Omega_{q}}{\partial\chi}=\mathrm{k_{0}\chi}\left(\sigma^{2}+\zeta^{2}+\delta^{2}\right)-k_{3}\left(\sigma^{2}-\delta^{2}\right)\zeta+\chi^{3}\left[1+{\rm{ln}}\left(\frac{\chi^{4}}{\chi_{0}^{4}}\right)\right]+(4k_{4}-d)\chi^{3}$
(38) $\displaystyle-$
$\displaystyle\frac{4}{3}d\chi^{3}{\rm{ln}}\Bigg{(}\bigg{(}\frac{\left(\sigma^{2}-\delta^{2}\right)\zeta}{\sigma_{0}^{2}\zeta_{0}}\bigg{)}\bigg{(}\frac{\chi}{\mathrm{\chi_{0}}}\bigg{)}^{3}\Bigg{)}+\frac{2\chi}{\chi_{0}^{2}}\left[m_{\pi}^{2}f_{\pi}\sigma+\left(\sqrt{2}m_{K}^{2}f_{K}-\frac{1}{\sqrt{2}}m_{\pi}^{2}f_{\pi}\right)\zeta\right]$
$\displaystyle-$
$\displaystyle\frac{\chi}{{\chi^{2}_{0}}}({m_{\omega}}^{2}\omega^{2}+{m_{\rho}}^{2}\rho^{2})=0,$
$\displaystyle\frac{\partial\Omega_{q}}{\partial\omega}=\frac{\chi^{2}}{\chi_{0}^{2}}m_{\omega}^{2}\omega+4g_{4}\omega^{3}+12g_{4}\omega\rho^{2}$
$\displaystyle-$ $\displaystyle\sum_{i=u,d}g_{\omega}^{i}\rho_{q,i}=0,$ (39)
$\displaystyle\frac{\partial\Omega_{q}}{\partial\rho}=\frac{\chi^{2}}{\chi_{0}^{2}}m_{\rho}^{2}\rho+4g_{4}\rho^{3}+12g_{4}\omega^{2}\rho$
$\displaystyle-$ $\displaystyle\sum_{i=u,d}g_{\rho}^{i}\rho_{q,i}=0,$ (40)
$\displaystyle\frac{\partial\Omega_{q}}{\partial\phi}=\frac{\chi^{2}}{\chi_{0}^{2}}m_{\phi}^{2}\phi+8g_{4}\phi^{3}$
$\displaystyle-$ $\displaystyle\sum_{i=s}g_{\phi}^{i}\rho_{q,i}=0,$ (41)
$\displaystyle\hskip
11.38092pt\frac{\partial\Omega_{q}}{\partial\Phi}=\bigg{[}\frac{-a(T)\bar{\Phi}}{2}-\frac{6b(T)(\bar{\Phi}-2{\Phi}^{2}+{\bar{\Phi}}^{2}\Phi)}{1-6\bar{\Phi}\Phi+4(\bar{\Phi}^{3}+\Phi^{3})-3(\bar{\Phi}\Phi)^{2}}\bigg{]}T^{4}-\sum_{i=u,d,s}\frac{2k_{B}TN_{C}}{(2\pi)^{3}}$
$\displaystyle\int_{0}^{\infty}d^{3}k\bigg{[}\frac{e_{q}\left(\frac{-(E_{i}^{*}(k)-{\nu_{i}}^{*})}{k_{B}T}\right)}{\left(1+e_{q}\left(\frac{-3(E_{i}^{*}(k)-{\nu_{i}}^{*})}{k_{B}T}\right)+3\Phi
e_{q}\left(\frac{-(E_{i}^{*}(k)-{\nu_{i}}^{*})}{k_{B}T}\right)+3\bar{\Phi}e_{q}\left(\frac{-2(E_{i}^{*}(k)-{\nu_{i}}^{*})}{k_{B}T}\right)\right)^{q}}$
$\displaystyle+\frac{e_{q}\left(\frac{-2(E_{i}^{*}(k)+{\nu_{i}}^{*})}{k_{B}T}\right)}{\left(1+e_{q}\left(\frac{-3(E_{i}^{*}(k)+{\nu_{i}}^{*})}{k_{B}T}\right)+3\bar{\Phi}e_{q}\left(\frac{-(E_{i}^{*}(k)+{\nu_{i}}^{*})}{k_{B}T}\right)+3\Phi
e_{q}\left(\frac{-2(E_{i}^{*}(k)+{\nu_{i}}^{*})}{k_{B}T}\right)\right)^{q}}\bigg{]}$
$\displaystyle=0,$ (42)

and

$\displaystyle\frac{\partial\Omega_{q}}{\partial\bar{\Phi}}=\bigg{[}\frac{-a(T)\Phi}{2}-\frac{6b(T)(\Phi-2{\bar{\Phi}}^{2}+{\Phi}^{2}\bar{\Phi})}{\mathrm{1-6\bar{\Phi}\Phi+4(\bar{\Phi}^{3}+\Phi^{3})-3(\bar{\Phi}\Phi)^{2}}}\bigg{]}T^{4}-\sum_{i=u,d,s}\frac{2k_{B}TN_{C}}{(2\pi)^{3}}$
$\displaystyle\int_{0}^{\infty}d^{3}k\bigg{[}\frac{e_{q}\left(\frac{-2(E_{i}^{*}(k)-{\nu_{i}}^{*})}{k_{B}T}\right)}{\left(1+e_{q}\left(\frac{-3(E_{i}^{*}(k)-{\nu_{i}}^{*})}{k_{B}T}\right)+3\Phi
e_{q}\left(\frac{-(E_{i}^{*}(k)-{\nu_{i}}^{*})}{k_{B}T}\right)+3\bar{\Phi}e_{q}\left(\frac{-2(E_{i}^{*}(k)-{\nu_{i}}^{*})}{k_{B}T}\right)\right)^{q}}$
$\displaystyle+\frac{e_{q}\left(\frac{-(E_{i}^{*}(k)+{\nu_{i}}^{*})}{k_{B}T}\right)}{\left(1+e_{q}\left(\frac{-3(E_{i}^{*}(k)+{\nu_{i}}^{*})}{k_{B}T}\right)+3\bar{\Phi}e_{q}\left(\frac{-(E_{i}^{*}(k)+{\nu_{i}}^{*})}{k_{B}T}\right)+3\Phi
e_{q}\left(\frac{-2(E_{i}^{*}(k)+{\nu_{i}}^{*})}{k_{B}T}\right)\right)^{q}}\bigg{]}$
$\displaystyle=0.$ (43)

## References

  * (1) D. Boyanovsky et al., Ann. Rev. Nucl. Part. Sci. 56, 441 (2006).
  * (2) M. Buballa et al., J. Phys. G 41, 123001 (2014).
  * (3) G. Baym et al., Rep. Prog. Phys. 81, 056902 (2018).
  * (4) B.-L. Li et al., Phys. Rev. D 99, 043001 (2019).
  * (5) T. Hinderer et al., Phys. Rev. D 81, 123016 (2010).
  * (6) F. Weber, Prog. Part. Nucl. Phys. 54, 193 (2005).
  * (7) R. Marty and J. Aichelin, Phys. Rev. C 87, 034912 (2013).
  * (8) K. Tsushima et al., Phys. Rev. C 59, 2824 (1999).
  * (9) K. Saito and A.W. Thomas, Phys. Lett. B 327, 9 (1994).
  * (10) B.J. Schaefer et al., Phys. Rev. D 76, 074023 (2007).
  * (11) T.K. Herbst et al., Phys. Lett. B 696, 58 (2011).
  * (12) D.P. Menezes et al., J. Phys. G 32, 1081 (2006).
  * (13) R. Stiele et al., Phys. Lett. B 729, 72 (2014).
  * (14) K. Fukushima, Phys. Lett. B 591, 277 (2004).
  * (15) C. Ratti, M.A. Thaler, W. Weise, Phys. Rev. D 73, 014019 (2006).
  * (16) M. Abu-Shady and H. M. Mansour, Phys. Rev. C 85, 5 (2012).
  * (17) R. Gatto and M. Ruggieri, Phys. Rev. D 78, 034015 (2011).
  * (18) B.J. Schaefer and M. Wagner, Prog. Part. Nucl. Phys. 62, 381 (2009).
  * (19) P. Kovacs, Zs. Szep, and Gy. Wolf, Phys. Rev. D 93 11, 114014 (2016).
  * (20) P. Kovacs, G. Kovacs, and F. Giacosa, Phys. Rev. D 106 11, 116016 (2022).
  * (21) P. Wang et al., Commun. Theor. Phys. 36, 71 (2001).
  * (22) P. Wang et al., Nucl. Phys. A 688, 791 (2001).
  * (23) M. Kumari and A. Kumar, Eur. Phys. J. Plus 136, 19 (2021).
  * (24) N. Chahal and A. Kumar, Chin. Phys. C 46, 6 (2022).
  * (25) C. Tsallis, J. Stat. Phys. 52, 479 (1998).
  * (26) Ya-Peng Zhao, Phys. Rev. D 101, 096006 (2020).
  * (27) J. Rozynek and G. Wilk, Eur. Phys. A 52, 13 (2016).
  * (28) T. Bhattacharyya et al., Eur. Phys. J. A 52, 30 (2016).
  * (29) Asmaa G. Shalaby, A. Phys. Pol. B 47, 5 (2016).
  * (30) G. Wilk and Z. Wlodarczyk, Eur. Phys. J. A 48, 161 (2012).
  * (31) B.-C. Li, Y.-Z. Wang, and F.-H. Liu, Phys. Lett. B 725, 352 (2013).
  * (32) L. Marques, J. Cleymans, and A. Deppman, Phys. Rev. D 91, 054025 (2015).
  * (33) B. De, Eur. Phys. J. A 50, 138 (2014).
  * (34) M. Rybczynski and Z. Wlodarczyk, Eur. Phys. J. C 74, 2785 (2014).
  * (35) G. Combe et al., Phys. Rev. Lett. 115, 238301 (2015).
  * (36) U. Tirnakli and E. P. Borges, Sci. Rep 6, 23644 (2016).
  * (37) L. J. L. Cirto, A. Rodriguez, F. D. Nobre et al., EPL 123, 30003 (2018).
  * (38) A. Adare et al. (PHENIX Collaboration), Phys. Rev. C 83, 064903 (2011).
  * (39) B. I. Abelev et al. (STAR Collaboration), Phys. Rev. C 75, 064901 (2007).
  * (40) G. Aad et al. (ATLAS Collaboration), New J. Phys. 13, 053033 (2011).
  * (41) K. E. A. Aamodt (ALICE Collaboration), Eur. Phys. J. C 71, 1655 (2011).
  * (42) V. E. A. Khachatryan (CMS Collaboration), J. High Energy Phys. 02, 041 (2010).
  * (43) I. Bediaga, E.M.F. Curado, and J.M. de Miranda, Physica A 286, 156 (2000).
  * (44) W.M. Alberico, A. Lavagno, and P. Quarati, Eur. Phys. J. C 12, 499 (2000).
  * (45) A. Lavagno, A.M. Scarfone, and P.N. Swamy, J. Phys. A 40, 8635 (2007).
  * (46) C. Beck, Eur. Phys. J. A 40, 267 (2009).
  * (47) G. Wilk and Z. Wlodarczyk, Phys. Rev. Lett. 84, 2770 (2000).
  * (48) A. Lavagno, P. Quarati, Phys. Lett. B 498, 47 (2001).
  * (49) T.S. Biro and G. Purcsel, Phys. Rev. Lett. 95, 162302 (2005).
  * (50) T.S. Biro, G. Purcsel, and K. Urmossy, Eur. Phys. J. A 40, 325 (2009).
  * (51) A. Lavagno, Phys. Lett. A 301, 13 (2002).
  * (52) W.M. Alberico et al., Physica A 387, 467 (2008).
  * (53) P. Quarati and A.M. Scarfone, Astrophys. J. 666, 1303 (2007).
  * (54) J.C. Carvalho et al., Astrophys. J. 696, L48 (2009).
  * (55) G. Livadiotis and D.J. McComas, Astrophys. J. 714, 971 (2010).
  * (56) M.P. Leubner, Astrophys. J. 632, L1 (2005).
  * (57) C. Tsallis, Introduction to nonextensive statistical mechanics: approaching a complex world (Springer Science, 2009).
  * (58) F. I. M. Pereira, R. Silva, and J. S. Alcaniz Phys. Rev. C 76, 015201 (2007).
  * (59) J. Roynek and G. Wilk 2009 J. Phys. G: Nucl. Part. Phys. 36 125108 (2009).
  * (60) P.H.G. Cardoso, T. Nunes da Silva, and A. Deppman, et al., Eur. Phys. J. A 53, 191 (2017).
  * (61) Shen et al., Adv. High Energy Phys. 2017 (2017), 4135329 (2017).
  * (62) Ya-Peng Zhao et al., Chinese Phys. C 45, 073105 (2021).
  * (63) Ya-Peng Zhao et al. Chinese Phys. C 47 053103 (2023).
  * (64) W.M. Alberico, A. Lavagno, and P. Quarati, Eur. Phys. J. C 12, 499 (2000).
  * (65) A. Drago, A. Lavagno, and P. Quarati, Physica A 344 472477 (2004).
  * (66) A. Lavagno et al, J. Phys. G: Nucl. Part. Phys. 37, 115102 (2010).
  * (67) A. Lavagno and D. Pigato, Eur. Phys. J. A 47, 52 (2011).
  * (68) E. Megias, D. P. Menezes, and A. Deppman, Physica A 421, 15 (2015).
  * (69) S.K. Tiwari et al. Eur. Phys. J. C 78, 938 (2018).
  * (70) J. Schechter, Phys. Rev. D 21, 3393 (1980).
  * (71) H. Gomm et al., Phys. Rev. D 33, 801 (1986).
  * (72) E.K. Heide, S. Rudaz, and P.J. Ellis, Nucl. Phys. A 571, 713 (1994).
  * (73) G. Carter, P.J. Ellis, and S. Rudaz, Nucl. Phys. A 603 , 367-386 (1996).
  * (74) P. Ko and S. Rudaz, Phys. Rev. D 50, 6877-6894 (1994).
  * (75) P. Papazoglou et al., Phys. Rev. C 59, 411 (1999).
  * (76) P. Wang et al., Phys. Rev. C 67, 015210 (2003).
  * (77) P. Wang et al., Nucl. Phys. A 705, 455 (2002).
  * (78) P. Wang et al., Nucl. Phys. A 744, 273 (2004).
  * (79) P. Costa et al., Symmetry 2, 1338 (2010).
  * (80) S. Roessner et al., Phys. Rev. D 75, 034007 (2007).
  * (81) S. Roessner et al., Nucl. Phys. A 814, 118 (2008).
  * (82) M. Fukugita, M. Okawa, and A. Ukava, Nucl. Phys. B 337, 181 (1990).
  * (83) D. H. E. Gross, Physica A 340, 76 (2004).
  * (84) D. H. E. Gross, Physica A 365, 138 (2006).
  * (85) B. C. Li, Y. Z. Wang, and F. H. Liu, Phys. Lett. B 725, 352 (2013).
  * (86) J. Cleymans et al., Phys. Lett. B 723, 351 (2013).
  * (87) M. D. Azmi and J. Cleymans, J. Phys. G 41, 065001 (2014).
  * (88) T. Kodama et al., 2005 Europhys. Lett. 70, 439 (2005).
  * (89) V. Garcia-Morales and J. Pellicer, Physica A  361, 161 (2006).
  * (90) A. M. Teweldeberhan, A. R. Plastino, H. C. Miller, Phys. Lett. A 343, 71 (2005).
  * (91) T. S. Biro and A. Jakovac, Phys. Rev. Lett. 94, 132302 (2005).
  * (92) A. Khuntia et al., Eur. Phys. J. A 52, 292 (2016).
  * (93) H. Mao et al., J. Phys. G: Nucl. Part. Phys. 37, 035001 (2010).