††thanks: Email<EMAIL_ADDRESS>Email<EMAIL_ADDRESS>
# Specific heat and its high-order moments in relativistic heavy-ion
collisions from a multiphase transport model

Ru-Xin Cao(曹汝鑫)  Shanghai Institute of Applied Physics, Chinese Academy of
Sciences, Shanghai 201800, China  Key Laboratory of Nuclear Physics and Ion-
beam Application (MOE), Institute of Modern Physics, Fudan University,
Shanghai 200433, China  School of Nuclear Sciences and Technology, University
of Chinese Academy of Sciences, Beijing 100049, China  Song Zhang(张松) ID Key
Laboratory of Nuclear Physics and Ion-beam Application (MOE), Institute of
Modern Physics, Fudan University, Shanghai 200433, China  Shanghai Research
Center for Theoretical Nuclear Physics， NSFC and Fudan University, Shanghai
200438, China  Yu-Gang Ma(马余刚) ID Key Laboratory of Nuclear Physics and Ion-
beam Application (MOE), Institute of Modern Physics, Fudan University,
Shanghai 200433, China  Shanghai Research Center for Theoretical Nuclear
Physics， NSFC and Fudan University, Shanghai 200438, China

###### Abstract

Energy dependence of specific heat extracted from temperature fluctuation of
Au + Au collisions at $\sqrt{s_{NN}}$ = 7.7 GeV to 200 GeV was investigated by
using a multiphase transport (AMPT) model. The results were compared with
those from other models and some differences at low $\sqrt{s_{NN}}$ were
found. To explain the above differences and describe the properties of the hot
dense matter at low $\sqrt{s_{NN}}$, a new quantity $C_{v}^{*}$ was derived
for describing specific heat in heavy-ion collisions. It was found that, by
using $C_{v}^{*}$ together with its high order moments (skewness and
kurtosis), thermal properties of the hot dense matter can be described and
different thermal properties with or without parton process can be clearly
distinguished. The proposed observable provides a way to learn the property of
QCD matter in heavy-ion collisions.

##  I Introduction

During the past few decades, amount efforts have been made on studies of the
hot dense quark matter created in relativistic heavy-ion collisions. Plenty of
evidences supported the existence of quark-gluon plasma (QGP) in relativistic
heavy-ion collisions and aroused the interest to explore the properties of hot
dense matter created under the extreme conditions of temperature and density
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. Specific heat was carried out as one of the
signals of phase transition in partonic level as well as nucleonic level and
helps to inspect thermal properties of nuclear matter. In partonic level, the
Lattice QCD predicted the phase transition as a crossover at zero baryon-
chemical potential $(\mu_{B}=0)$ but while a first order phase transition
could occur at finite baryon density [11, 12, 13, 14, 15, 16, 17]. In
nucleonic level, a liquid gas phase transition can occur at sub-saturation
density and moderate temperature [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29].

According to phase transition theory, the long-range correlation diverges
rapidly when a thermodynamic system evolves close to critical point. Here it
is relevant to the degree of freedom of the quarks and gluons, rather than
mesons and baryons for the matter created in the collisions. As a response to
the system perturbation the specific heat could diverge when the system
evolves close to the critical point. In statistical physics the heat capacity
can be associated with the fluctuation of temperature. And many works tried to
extract heat capacity from the event-by-event distribution of ensemble
temperature [16, 30, 31, 32, 33, 34, 35, 12, 13, 14, 15, 11, 23, 25]. At near
the critical point, the specific heat $C_{V}$ can be expressed as a power law
function of temperature deviated from critical temperature ($T_{c}$), i.e.
$C_{V}\propto|T-T_{c}|^{-\alpha}$ with $\alpha$ as a critical exponent.
Furthermore, although energy dependence of the specific heat $C_{V}$ has been
studied in some theoretical works [33, 34, 36], event-by-event fluctuation of
the specific heat itself should be also paid more attention for investigating
its statistical properties as a potential probe to the QCD critical point.

To this end, a multiphase transport (AMPT) model [37, 38, 39] was adopted to
simulate Au + Au collisions at various energies of $\sqrt{s_{NN}}$ = 7.7,
11.5, 14.5, 19.6, 27, 39, 62.4, and 200 GeV. Energy dependence of $C_{V}$ was
obtained and compared with results extracted through the same temperature
fluctuation method by other models, namely the Hadron Matter (HM) [33], Hadron
Matter via Quark Gluon Matter [33], Quark Gluon Matter (QGM) [33], Hadron
Resonance Gas (HRG) models as well as the data deduced from the STAR
Collaboration [34, 36]. However, it was found that the specific heat obtained
from the AMPT model shows obvious depression at low $\sqrt{s_{NN}}$, similar
to some other models [34, 36, 33]. This insensitivity phenomenon of
temperature fluctuation in the AMPT model drove our present study to propose
another derivation of heat capacity from the basic definition with some
assumptions. The derivation in this study gives an effective specific heat,
denoted as $C_{v}^{*}$, which is expressed by the characteristic event’s
kinetic quantities, and plots the event-by-event distribution to obtain the
skewness and kurtosis to describe the statistical properties of $C_{v}^{*}$.
The mean value of new derived $C_{v}^{*}$ together with its skewness and
kurtosis give insight of the evolution of thermodynamical properties of the
emitted particles and describe energy dependence of specific heat, which could
be performed in experiments.

The paper is organized as follows: in Sec. II, a brief introduction of a
multiphase transport model was given, then the results on the energy
dependence of specific heat from the event-by-event temperature fluctuation in
Au + Au collisions simulated by the AMPT were present, and the previous
results from other models as well as the STAR data were compared. By the above
comparison, the difference and the uncertainty of specific heat at low
$\sqrt{s_{NN}}$ in our AMPT calculation were analyzed. In Sec. III, a brief
qualitative explanation for the differences was presented, and a formula
expressed by kinetic quantities called $C_{v}^{{}^{*}}$ from the definition of
specific heat was proposed. Sec. IV gave the results for mean value, skewness,
kurtosis of event-by-event distributions for $C_{v}^{*}$, from which the
meaning of those moments was explained by emitted particles’ thermodynamic
properties. Finally a brief summary was given in Sec. V.

##  II AMPT results and comparison

###  II.1 Introduction to AMPT

A multiphase transport (AMPT) model [37, 38, 39] is composed of four stages to
simulate relativistic heavy-ion collisions. It has successfully described
various phenomena at RHIC and LHC energies and becomes a well-known event
generator. The AMPT has two versions: String Melting (SM) and Default. In SM
version, Heavy Ion Jet Interaction Generator (HIJING) [40, 41] is used to
simulate the initial conditions, then Zhang’s Parton Cascade (ZPC) [42] is
used to describe interactions for partons which are from all of hadrons in the
HIJING but spectators, after which a simple Quark Coalescence Model describes
hadronization process, finally A Relativistic Transport (ART) model [43]
simulates hadron re-scattering process. The Default version of AMPT only
conducts the minijet partons in partonic scatterings via ZPC and uses the Lund
string fragmentation to perform hadronization.

AMPT model [37, 39] can describe the $p_{T}$ spectrum and energy dependence of
identified particles such as pion, kaon, $\phi$, proton and $\Omega$ produced
in heavy-ion collisions [38, 44, 45], as well as the collective flows and
temperature during evolution etc [46, 47, 34, 48]. Chiral and magnetic related
anomalous phenomena can also be described by the AMPT model [49, 50, 51, 52,
53, 54]. More details of the model description and selection of the set for
parameters can be found in Refs. [37, 38, 39]. The present study used both
versions to simulate Au + Au collisions at $\sqrt{s_{NN}}$ = 7.7, 11.5, 14.5,
19.6, 27, 39, 62.4, and 200 GeV in the impact parameter range of 0 - 4.7 fm
which is corresponding to centrality 0 - 10%.

###  II.2 $C_{V}$ from temperature fluctuation

Firstly, an event-by-event fluctuation method was used to study the specific
heat. As noted in Refs. [34, 33], since the lower multiplicity of particles
caused relatively large error in fitting $p_{T}$ spectrum, we then merged some
events into one new event, denoted as “linked event”. The combined events’
number is determined by ensuring the multiplicity in one “linked event” more
than 1000, thus we increased mean multiplicity in one “event” more than 1000,
as other researchers did in Ref. [34]. And as shown in Fig. 8 of Ref. [34], by
comparing the results in “linked events” and in “random mixed events” in which
each particle from different events (plotted by blue and red line), we can see
“linking events” from different original event sets only modify fluctuation
results slightly.

The effective temperature, $T_{eff}$, of a $\pi^{+}$ system at final state is
obtained via fitting $p_{T}$ spectrum by using the exponential distribution
[38, 55, 56]:

$\displaystyle\begin{aligned}
\frac{1}{p_{T}}\frac{dN}{dp_{T}}=Ae^{\frac{-p_{T}}{T_{eff}}}.\end{aligned}$
(1)

$\pi^{+}$ system is chosen for the calculations with kinetic window, $p_{T}<2$
GeV/$c$ and $-1<y<1$. $T_{eff}$ stands for the slope of $p_{T}$ spectrum and
consists of two parts， i.e. the kinetic temperature $T_{kin}$ and contribution
from radial flow $\langle\beta_{T}\rangle$ [34],

$\displaystyle T_{eff}=T_{kin}+f(\langle\beta_{T}\rangle),$ (2)

where the kinetic freeze-out temperature, $T_{kin}$, characterizes thermal
motion of emitted particles, and $f(\langle\beta_{T}\rangle)$ reflects
transverse radial flow contribution. As we are discussing $\pi^{+}$ system,
the relation in Eq. (2) can be approximately written as [34],

$\displaystyle T_{eff}\approx T_{kin}+m_{0}\langle\beta_{T}\rangle^{2},$ (3)

where $m_{0}$ is the mass of $\pi^{+}$.

For a system in equilibrium, its heat capacity ($C$) is related to the
fluctuation of event-by-event temperature distribution [35, 14, 11, 15],

$\displaystyle\begin{aligned} P(T)\propto e^{-\frac{C}{2}\frac{(\Delta
T)^{2}}{\langle T\rangle^{2}}},\end{aligned}$ (4)

where $\langle T\rangle$ is the mean temperature value and $\Delta T=T-\langle
T\rangle$ is the deviation of temperature from its mean value. The expression
can be further derived as [34]:

$\displaystyle\begin{aligned} \frac{1}{C}=\frac{\langle
T_{kin}^{2}\rangle-\langle T_{kin}\rangle^{2}}{\langle
T_{kin}\rangle^{2}}\approx\frac{\langle T_{eff}^{2}\rangle-\langle
T_{eff}\rangle^{2}}{\langle T_{kin}\rangle^{2}}.\end{aligned}$ (5)

When multiplicity $N$ is taken into account, we can get specific heat per
particle ($C_{V}$) as the following

$\displaystyle\begin{aligned} C_{V}=\frac{C}{N}.\end{aligned}$ (6)

In this section we choose $T_{kin}$ to extract $C_{V}$ since the variances of
$T_{kin}$ and $T_{eff}$ are close enough to each other as mentioned in Ref.
[34]. Figure 1 shows the $p_{T}$ spectra in “linked events” fitted by Eq. (1),
from which the effective temperature was obtained as shown in Fig. 2(a).
Meanwhile, in the present work we calculated transverse flow velocity
$\beta_{T}$ by
$\beta_{T}=\frac{\bm{\beta}_{\bm{p_{T}}}\cdot\bm{\rho}}{|\bm{\rho}|}$, where
$\bm{\beta}_{\bm{p_{T}}}$ is the transverse velocity and $\bm{\rho}$ is the
transverse coordinates of particles at kinetic freeze-out state. In Refs. [34,
57] where a Blast-Wave Model was used to obtain a fit parameter also denoted
by $\beta_{T}$. The two definitions give the same dimension of $\beta_{T}$ and
are both on a super-surface where particles reach kinetic freeze-out status.
So the two $\beta_{T}$ are conceptually consistent. As mentioned in Ref. [34],
fluctuation of $\beta_{T}$ could be dominant for small, asymmetric and non-
central collisions. However, our work focus on Au + Au system, the fluctuation
from $\beta_{T}$ is almost negligible.

As we mentioned before, we used the Boltzmann distribution to extract
effective temperature from $p_{T}$ spectra as shown in Fig. 1. However, it is
noticed that as $p_{T}$ increases, the Boltzmann distribution becomes poor to
fit the spectra. To ensure our fitting results comparable and verifiable to
results in Refs. [34, 33], however, we limit the fitting range of $p_{T}$
below 1 GeV/$c$ with the Boltzmann distribution, which was represented by the
lines in Fig. 1. Of course, some checks have been done. For examples, one way
is extending the $p_{T}$ fit range to 3 GeV/c and another way is just checking
the non-thermal tail of high $p_{T}$ from 1 GeV/c to 3 GeV/c with the same
Boltzmann distribution. From the above fit procedures, we can roughly get the
lower ($p_{T}<$ 1 GeV/c) and upper (1 GeV/c $<p_{T}<$ 3 GeV/c) limits of the
event-by-event temperature distribution. Further we can get the lower and
upper bounds for heat capacity. It was found that the different $p_{T}$ range
does not change the energy dependence of $C_{V}$ but do induce an uncertainty
of about 30% of the $C_{V}$ values.

On the other hand, there exists a few alternative models to fit $p_{T}$
spectra, eg. with the Tsallis-Blast-Wave model to describe $p_{T}$ spectra in
a wide $p_{T}$ range up to 3 GeV/$c$ shown in Ref. [58, 59]. The Tsallis
distribution introduced a new parameter $q$ to represent how the system
deviates from the equilibrium state. And it can be seen in Ref. [58] that in
our range for $p_{T}<$1 GeV/$c$, $q-1$ is very close to zero. The fact that
$q-1$ tends to $0$ means the Tsallis distribution degrades into the Boltzmann
distribution, indicating that our fits using the Boltzmann distribution is
reasonable. Further, in Ref. [60], the Tsallis distribution was used to
investigate thermal parameters in small system ($p+p$ at $\sqrt{s_{NN}}=7$
TeV) and fitted well in the PYTHIA frame. Of course, adopting a new Tsallis-
Blast-Wave fit to deduce the heat capacity is naturally an interesting topic,
which deserves further investigation in future.

Figure 1: The $p_{T}$ spectrum for the merged events with $-1<y<1$ at
$\sqrt{s_{NN}}=7.7-200$ GeV in $0-10\%$ central Au + Au collisions by using
the String Melting AMPT. The line represents the Boltzmann fit to the $p_{T}$
spectra below 1 GeV/c.

After removing transverse radial flow contribution from Eq. (3), the
distribution of event-by-event $T_{kin}$ is presented in Fig. 2(b) at the
centrality of $0-10\%$ and $\sqrt{s_{NN}}$ = 7.7 - 200 GeV. It is clearly seen
that $T_{kin}$ is much lower than $T_{eff}$ because the transverse radial flow
is taken off. Via fitting the $T_{kin}$ distribution by Eq. (5) the heat
capacity per particle $C_{V}$ can be extracted.

Figure 2: The event-by-event distributions of effective temperature (a) and
kinetic temperature (b) at $\sqrt{s_{NN}}=7.7-200$ GeV in $0-10\%$ central Au
+ Au collisions by using the String Melting AMPT model, the lines are fitted
by Gaussian distribution.

Fig. 3 shows different model results together with the STAR data by
temperature fluctuation method. The extracted $C_{V}$ from the STAR data
decreases with the increasing of $\sqrt{s_{NN}}$ [34]， and it shows a
significant rising trend when $\sqrt{s_{NN}}$ goes down to $40$ GeV and below,
which indicates that the system might enter critical region. The HRG model
[34] shows similar energy dependence of $C_{V}$ to the STAR Au + Au results.
But other results from the HM [33], the QGM [33], the AMPT [34], as well as
this work employing AMPT with removing the radial flow effect give an almost
constant value of $C_{V}$ in a wide range of $\sqrt{s_{NN}}$. These results
demonstrated a significant model dependence at low $\sqrt{s_{NN}}$ but are
consistent with each other at high $\sqrt{s_{NN}}$ [34, 36, 33].

To understand the difference between the AMPT simulation results and the
enhancement structure from the STAR data or the HRG model [34] at low
$\sqrt{s_{NN}}$, a few issues may be related as indicated in Refs. [34, 33]:

(1) The finite particle multiplicity, choice of kinetic window of $p_{T}$ and
fit parameters in Eq. (1) provide uncertainties of $C_{v}$, which was
mentioned in Refs. [34, 33] as well. To check our fitting results on $p_{T}$
spectrum and ensure our fitting results comparable to those in Ref. [34], we
chose the same $p_{T}$ window as Ref. [34].

(2) In many theoretical derivations, the system evolves in perfect
thermodynamic conditions, which means that the system volume is fixed and
multiplicity should be high enough. However, these conditions can hardly be
satisfied in real heavy-ion collision experiments, as noted in Ref. [33];

Two reasons mentioned above are both independent of models. Besides, the model
dependence should also be commented, i.e., different model frames provide
respective additional uncertainty, which could be caused by different
mechanisms used in models, even the results are calculated in the same kinetic
window of transverse momentum and rapidity.

In the following, we try to discuss the insensitivity of temperature
fluctuation in the present AMPT framework.

Figure 3: The energy dependence of specific heat $C_{v}$ = $\frac{C}{N}$ is
calculated via the temperature fluctuation method: the AMPT model with
removing radial flow effect (this work), the results using the STAR 0-5%
centrality Au-Au data (Basu et al. 2016 [34]), the HRG model (Basu et al. 2016
[34]), the AMPT model with flow included (Basu et al. 2016 [34]), the HM model
(Sa et al. 2007 [33]), the HM via the QGM model (Sa et al. 2007 [33]), and the
QGM model (Sa et al. 2007 [33]) .

##  III Derivation of a new formula of heat capacity: $C_{v}^{*}$

To avoid the model dependence of extracted specific heat, the definition of
heat capacity according to Ref. [61] can be written as

$\displaystyle\begin{aligned} C=\left(\frac{\partial E_{th}}{\partial
T}\right)_{V},\end{aligned}$ (7)

where $E_{th}$ is the energy of thermal motion in a researched system and $T$
is the temperature with a fixed volume of the system. Actually the definition
of the heat capacity by Eq. (7) and Eq. (5) are equivalent. In Ref. [61] there
is nothing different but the one in Ref. 7 is defined from differential of
enthalpy in Eq. (14.6) of Chapter 15 and the one in Ref. 5 is derived from
fluctuation in Eq. (111.6) of Chapter 12.

The definition in Eq. (7) includes following assumptions:

(1) The energy of $\pi^{+}$ system at final state – $E$ – has positive
correlation to $\sqrt{s_{NN}}$, which means for each collision at final state,
total energy of $\pi^{+}$ system evolves along a certain curve with
$\sqrt{s_{NN}}$ [62];

(2) Estimated phase volume can be approximately calculated as charged particle
multiplicity $N$, which means that we can calculate specific heat per particle
$C_{v}=\frac{C}{N}$ and the merged events are ensured to have enough
multiplicity (here over one thousand) [34, 33];

(3) The $\pi^{+}$ system evolves to the kinetic freeze-out stage after the
hadronic interaction. That means the the energy contributed by system’s
internal interaction (mainly the hadronic interaction) can be ignored when we
calculated the heat capacity.

(4) The total energy of system is mainly dominated by transverse momentum with
$\pi^{+}$ at mid-rapidity, which can be concluded from the transverse velocity
distribution for $\pi^{+}$ as shown in Fig. 4, most $\pi^{+}$ are distributed
in the region over $\beta_{p_{T}}>0.8c$. So $\pi^{+}$’s energy can be written
as a function of $p_{T}$. Considering the radial flow, the total energy can be
written as $E=E_{th}+E_{R}$ where $E_{th}$ is internal energy from thermal
motion and $E_{R}$ is energy from collective radial flow [34]. In
thermodynamical viewpoint, the heat capacity should only be correlated to
$E_{th}$. We can take a quick glance for the radial collective motion and
thermal motion of the system via giving distributions of radial flow
$\beta_{T}$ and thermal velocity $\beta_{th}$ =
$\frac{\bm{\beta}_{\bm{p_{T}}}\cdot\bm{\rho_{T}}}{|\bm{\rho_{T}}|}$ where
$\bm{\rho_{T}}$ and $\bm{\rho}$ constitute a right-handed orthogonal system.
As shown in Fig. 4, the pions are moving thermally in a expanding fireball
with a collective radial velocity. To avoid complex differential term of
$E_{R}$, we firstly use total energy $E$ to derive a formula of an effective
heat capacity. The following calculation will show us that, though the derived
effective heat capacity includes contribution from $E_{R}$, the radial flow
still can not dominate the energy dependence and fluctuation properties of
heat capacity. This fact ensures our calculation and discussion on new derived
heat capacity reasonable.

Figure 4: The $\beta_{p_{T}}$, $\beta_{T}$ and $\beta_{th}$ distribution for
$\pi^{+}$ system at $\sqrt{s_{NN}}$ = 7.7 and 200 GeV in 0-10% central Au + Au
collisions by using the String Melting AMPT model.

The total energy of a $\pi^{+}$ system can be written as [37, 46]

$\displaystyle\begin{aligned}
E=\sum_{i=1}^{N}E_{i}=\sum_{i=1}^{N}\left(p_{Ti}^{2}+p_{zi}^{2}+m_{0}^{2}\right)^{\frac{1}{2}}.\end{aligned}$
(8)

By calculating the variance and correlation coefficient of $T_{kin}$ and
radial flow, we can prove that in the following derivation of $C_{v}$, the
$E_{th}$ can be replaced by total energy $E$ in Eq. (8) reasonably.

Here we discuss the contribution in energy from radial flow. It should be
noted that we did not remove $f(\langle\beta_{T}\rangle)$ from $E$ directly
because our $\beta_{T}$ is not exactly same as the one in the Blast-Wave Model
in energy composition (though in temperature calculation our $\beta_{T}$
approximately equals to $\beta_{T}$ from the Blast-Wave model [55, 57]). Thus
we can not simply write the energy contributed by radial flow as form of non-
relativistic kinetic energy [57, 63, 64].

In Eq. (3), the variance of $T_{eff}$ can be written as,

$\displaystyle\begin{aligned}
\sigma^{2}_{T_{eff}}\approx\sigma^{2}_{T_{kin}}+m_{0}^{2}\sigma^{2}_{\langle\beta_{T}\rangle^{2}}+2m_{0}\mathrm{Cov}\left(T_{kin},\langle\beta_{T}\rangle^{2}\right).\end{aligned}$
(9)

Fig. 5 shows the variance of $T_{eff}$ and each term in the right side of Eq.
(9). We see that the variance from $T_{kin}$ dominates the one for $T_{eff}$
while the one from $\langle\beta_{T}\rangle^{2}$ is less with one order of
magnitude, thus the variance of $\langle\beta_{T}\rangle$ could be negligible
for properties of temperature. In order to clarify that the radial flow term
is independent from temperature term in our derivation and calculation, we
investigate the linear correlation coefficient
$\eta=\mathrm{Cov}(T_{kin},\langle\beta_{T}\rangle^{2})/(\sigma_{T_{kin}}\sigma_{\langle\beta_{T}\rangle^{2}})$
which is also presented in Fig. 5, here $\eta\leq 0.2$ means that they are
linear independent. In addition the derivative of $T_{eff}$ over $T_{kin}$ is
approximately a constant, then we can use $T_{eff}$ instead of $T_{kin}$ for
differential term from radial flow.

Figure 5: The variance of each term in temperature from the String Melting
AMPT model, and the correlation coefficient between $T_{kin}$ and
$\beta_{T}^{2}$ in different versions.

According to our assumptions, a $\pi^{+}$ system evolves along a continuous
path on $\sqrt{S_{NN}}$. Each system with its evolution path corresponds to
their own kinetic parameters like $E$, $T$ and $\langle p_{T}\rangle$, when
$\sqrt{s_{NN}}$ changes, the total energy $E$ of $\pi^{+}$ system changes to
$E+\Delta E$ on the continuous path. Similarly, when the $p_{T}$ spectrum
changes, resulting $p_{T}$ of each $\pi^{+}$ in one system changes to
$p_{T}+\Delta{p_{T}}$. Here we see the advantage of the merged events, which
ensured similar multiplicity when $E$ changes. That means for each $\pi^{+}$
on the continuous evolution path we can find a closest new $\pi^{+}$ in phase
space with $p_{T}+\Delta{p_{T}}$. In this sense, $E$ and $p_{T}$ for each
$\pi^{+}$ can be seen as a continuous function of system temperature.

Under the above four assumptions, especially considering the fact that the
energy only depends on $p_{T}$ in mid-rapidity, by using Eq. (7) and Eq. (8)
we get

$\displaystyle\begin{aligned} \frac{\partial E}{\partial
T}|_{E\left(T\right)}&=\sum_{i}^{N}\frac{\partial}{\partial
T}\left(p_{Ti}^{2}+p_{zi}^{2}+m_{0}^{2}\right)^{\frac{1}{2}}\\\
&=\sum_{i}^{N}\frac{\partial E_{i}}{\partial p_{Ti}}\frac{\partial
p_{Ti}}{\partial T}=\sum_{i}^{N}\frac{1}{2}\frac{2p_{Ti}}{E_{i}}\frac{\partial
p_{Ti}}{\partial T}\\\ &=\sum_{i}^{N}\left(\beta_{p_{Ti}}\frac{\partial
p_{Ti}}{\partial T}\right)|_{E\left(T\right)}.\end{aligned}$ (10)

The sum can be converted to the integral pattern with help of Riemann Sum. Let
$f\left(x_{i}\right)=F\left(x_{i}\right)/{\Delta x_{i}}$, for equidistant
division, the Riemann integral can be written as

$\displaystyle\begin{aligned}
\int_{a}^{b}f\left(x\right)g\left(x\right)dx\approx\sum_{i}^{\left(b-a\right)/{\Delta
x}}\frac{F\left(x_{i}\right)}{\Delta x}g\left(x_{i}\right)\Delta
x.\end{aligned}$ (11)

It should be noted here $N_{x}=\frac{\left(b-a\right)}{\Delta x}$ stands for
the number of division in Riemann integral which is different from
multiplicity $N_{multi}$, but we can approximately let $N_{x}=N_{multi}$, so
that the definition of $C_{v}$ expressed by sum can be related to integral
over temperature. As mentioned in Ref. [33] and Ref. [34], merged events
ensured enough high multiplicity and then small error for fitting.

Using Eq. (11), Eq. (10) can be expressed as follows,

$\displaystyle\begin{aligned}
&\sum_{i}^{N_{multi}}\beta_{p_{Ti}}\frac{\partial p_{Ti}}{\partial
T}|_{T}=\sum_{i}^{\left(b-a\right)/{\Delta x}}\frac{\beta_{p_{Ti}}}{\Delta
T}\frac{\partial p_{Ti}}{\partial T}\Delta T\\\
&\approx\int_{a}^{b}\frac{\beta_{p_{Ti}}}{dT}\frac{\partial p_{Ti}}{\partial
T}dT=\left.\left(\frac{\beta_{p_{Ti}}}{dT}\frac{\partial p_{Ti}}{\partial
T}\right)\right|_{T={T_{0}}}\left(b-a\right)\\\
&\left.\approx\tilde{\beta_{p_{T}}}\frac{\partial{\tilde{p_{Ti}}}}{\partial
T}\right|_{T={T_{0}}}\frac{\left(b-a\right)}{\Delta
T}=N_{multi}\left.\left(\tilde{\beta_{p_{T}}}\frac{\partial{\tilde{p_{Ti}}}}{\partial
T}\right)\right|_{T={T_{0}}}.\end{aligned}$ (12)

Here $T_{0}$ means the temperature in characteristic event
${}^{\prime}\xi^{\prime}$, corresponding to the ${}^{\prime}\xi^{\prime}$
noted in Mean Value Theorem of Integrals.

The Mean Value Theorem of Integrals says: Let
$f:\left[a,b\right]\rightarrow\mathbb{R}$ be a continuous function,
$g:\left[a,b\right]\rightarrow\mathbb{R}$ be integrable and constant-sign,
then there exists $\xi\in\left(a,b\right)$ so that
$\int_{a}^{b}f\left(x\right)g\left(x\right)dx=f\left(\xi\right)\int_{a}^{b}g\left(x\right)dx$.

According to Mean Value Theorem of Integrals , we can prove the existence of
$\tilde{\beta_{p_{T}}}$ and $\tilde{p_{Ti}}$, these two parameters reveal the
evolution property and energy dependence of system. Hence the derivation gave
the specific heat expressed by kinetic quantities:

$\displaystyle\begin{aligned}
C_{v}^{T_{0}}=\frac{C}{N_{multi}}\approx\tilde{\beta_{p_{T}}}\left({T_{0}}\right)\left.\frac{\partial\tilde{p_{T}}}{\partial
T}\right|_{T_{0}}.\end{aligned}$ (13)

In Eq. (13) and Eq. (10), the heat capacity for ensemble can be expressed as a
combination of kinetic measurements in a characteristic event whose
$T={T_{0}}$.

The above equation gave a formula for calculating specific heat approximately
by transverse velocity, transverse momentum and temperature in a unique event
${}^{\prime}\xi^{\prime}$, i.e. by a derivative of transverse momentum to
temperature in a unique event ${}^{\prime}\xi^{\prime}$. It should be noted
here that $C_{v}^{T_{0}}$ apparently includes radial flow term by $\beta_{T}$,
but $\beta_{T}$ decreases at low $\sqrt{s_{NN}}$ [55] and
$f(\langle\beta_{T}\rangle)$ is smaller than $T_{kin}$, the dominant
contribution is actually from thermal motion. That is why in Eq. (7) we can
use $E$ instead of $E_{th}$.

However, the first mean value theorem for integrals can only provide the
existence of ${}^{\prime}T_{0}^{\prime}$, it can not find the accurate value
to calculate the $C_{v}^{T_{0}}$. Therefore the present study needs to further
use the mean value of kinetic quantities instead of unique one to give
accessible results that can be measured and computed from data:

$\displaystyle\begin{aligned}
{C_{v}^{*}}=\overline{\beta_{p_{T}}}\left.\frac{\partial\overline{p_{T}}}{\partial
T}\right|_{T=\overline{T_{eff}}}.\end{aligned}$ (14)

We should emphasize that the definition in Eq. (14) is a compromise for the
real calculation of $C_{v}^{T_{0}}$. Since $\xi$ in the Mean Value Theorem for
Integrals is unsolvable, we use the mean value of $\beta_{p_{T}}$ and
$T_{eff}$ instead of values in single event to make sure $C_{v}^{T_{0}}$ can
be approximately calculated in practice.

According to the definition in Eq. (7), we should use $T_{kin}$ as $T$ in
calculation of $\overline{p_{T}}$ in Eq. (14). However, this formula involving
radial flow can not give identical results of
$\frac{\partial\overline{p_{T}}}{\partial T_{kin}}$. Fortunately, as
demonstrated in Fig. 5, it can be calculated by the AMPT result that
derivatives for $T_{kin}$ and $T_{eff}$ have strongly linear dependence, we
can directly use $T_{eff}$ instead of $T_{kin}$ in
$\frac{\partial\overline{p_{T}}}{\partial T}$. Now we can conveniently
calculate $\frac{\partial\overline{p_{T}}}{\partial
T}|_{T=\overline{T_{eff}}}$ based on $p_{T}$ spectrum [55, 65, 66, 67],

$\displaystyle\begin{aligned}
&\overline{p_{T}}=\frac{\int_{p_{min}}^{p_{max}}p_{T}^{2}F\left(p_{T}\right)dp_{T}}{\int_{p_{min}}^{p_{max}}p_{T}F\left(p_{T}\right)dp_{T}}=2T_{eff}+\\\
&\frac{p_{min}^{2}e^{-p_{min}/T_{eff}}-p_{max}^{2}e^{-p_{max}/T_{eff}}}{\left(p_{min}+T_{eff}\right)e^{-p_{min}/T_{eff}}-\left(p_{max}+T_{eff}\right)e^{-p_{max}/T_{eff}}}.\end{aligned}$
(15)

Considering that the mean value and unique value are obviously different,
which directly connected to our explanation on the statistical parameters, we
will discuss the topic in the next section.

##  IV Results of $C_{v}^{*}$ and its high order moments

Using the derived formula (14) and (LABEL:eq:<p_T>), $C_{v}^{*}$ can be
calculated with data of $\pi^{+}$ system in Au + Au collisions at different
$\sqrt{s_{NN}}$ simulated by the AMPT. The event-by-event distribution of
$C_{v}^{*}$ presents the Gaussian-like distribution as shown in Fig. 6, from
which the mean values ($\langle C_{v}^{*}\rangle$) at each $\sqrt{s_{NN}}$ can
be obtained.

Figure 6: (a) The event-by-event $C_{v}^{*}$ distribution for $\pi^{+}$ system
in $0-10\%$ central Au+Au collisions at $\sqrt{s_{NN}}$ = 7.7, 11.5, 14.5,
19.6, 27, 39, 62.4, and 200 GeV, respectively, by using the String Melting
AMPT. The lines are fitted by Gaussian distribution; (b) The extracted
$\langle C_{v}^{*}\rangle$ from the String Melting AMPT model calculation in
red line, compared to the $C_{v}$ from $T_{kin}$ fluctuation which is
represented by the black line, corresponding to the black star marker in Fig.
3. Note that the $C_{v}$ value from $T_{kin}$ fluctuation are scaled by twice.

Figure 6(b) shows $\left\langle C_{v}^{*}\right\rangle$ from the AMPT
calculations, which displays a slightly decreasing trend with the increasing
of $\sqrt{s_{NN}}$. Compared to the results in Ref. [33], the energy
dependence of $\left\langle C_{v}^{*}\right\rangle$ is similar to the results
of QGM, HM and HM via QGM in Fig. 3, even though the values of $C_{V}$ and
$\left\langle C_{v}^{*}\right\rangle$ are slightly different. In Ref. [33],
the specific heat from QGM is much different from results in pQCD
thermodynamic method and approximate pure gauge theory, which should be
studied further to confirm the discrepancy contributed by classical and
quantum statistics, the latter two theoretical calculations both gave larger
$C_{V}(10\sim 15)$ at low $\sqrt{s_{NN}}$ than other models shown in Fig. 3.
So they suggested using $p_{T}$ distribution of $\pi^{+}$ in RHIC single
events to extract $C_{V}$ according to the same statistical method as Ref.
[33]. For our calculation of $C_{v}^{*}$, we suggest extracting transverse
velocity and effective temperature from $p_{T}$ spectrum. But as far as the
uncertainty for present experimental data, if one considers to extract
parameters for $C_{v}^{*}$ which we discussed above, we suggest to improve
experimental precision to ensure smaller uncertainty for fitting parameters.
And if possible, the direct particle identification, transverse velocity
measurement and $p_{T}$ spectra for merged events of charged $\pi$ system in
experiments are also recommended, so that we can check our methods and fit
experimental results as accurate as possible.

To clearly illustrate the sensitivity of $C_{v}^{*}$ to $\beta_{p_{T}}$ and
$T_{eff}$, we plot $C_{v}^{*}$ as functions of $\beta_{p_{T}}$ and $T_{eff}$
by using Eqs. (14) and (LABEL:eq:<p_T>) in Fig. 7. An inflection point exists
at around $T_{eff}=100$ MeV, which depends on the fitting limits of $p_{T}$
spectrum. The combination of linear dependence of $\beta_{p_{T}}$ and non-
monotonic dependence of $T_{eff}$ results in that $C_{v}^{*}$ is sensitive to
the range of the fitting parameters. Along an identified curve of $C_{v}^{*}$
in Fig. 7, a shifting fitting result may cause different trend of energy
dependence of $C_{v}^{*}$, which has to be studied further, especially near
the QCD critical point.

Figure 7: Theoretical $C_{v}^{*}$ distribution as a two-dimensional function
of $\beta_{p_{T}}$ and $T_{eff}$.

In Fig. 6(b) the mean values of $C_{v}^{*}$ are smaller than those extracted
from the HRG model, HM model as well as QGM model [33]. In the first mean
value theorem for integrals, mean value deviates the ideal $C_{v}^{T_{0}}$.
That means only if one can accurately find the evolution curve of a
characteristic event corresponding to temperature $T_{0}$, the real
$C_{v}^{T_{0}}$ can agree with results extracted by temperature fluctuation
from the STAR’s data especially at low $\sqrt{s_{NN}}$ [34]. It is interesting
to see that $\langle C_{v}^{*}\rangle$ remains the rising trend with the
decreasing of $\sqrt{s_{NN}}$, which suggests that in the event-by-event
distribution of $C_{v}^{*}$ we can still find special properties similar to
sharp decreasing trend with increasing of colliding energy [47].

As a comparison, the behaviour of $C_{v}^{*}$ with parton process or only with
hadron process are different. Though the rising trend of average $C_{v}^{*}$
is depressed, the events with parton process provides higher $C_{v}^{*}$
similar to the cases of pQCD and other models [34, 33]. That means that more
events with higher $C_{v}^{*}$ emerge at low $\sqrt{s_{NN}}$ in the SM-AMPT
case. Further, we can investigate the fluctuation behaviour of $C_{v}^{*}$ by
event-by-event distribution as we did in temperature distribution. And indeed
we can see nontrivial fluctuation behaviour of $C_{v}^{*}$ when taking parton
interaction into account.

Based on the event-by-event distribution of $C_{v}^{*}$, skewness and kurtosis
of $C_{v}^{*}$ can be obtained in the AMPT model with the String Melting and
Default version, respectively. The skewness and kurtosis are defined as
$\frac{\mu_{3}}{\sigma^{3}}$ and $\frac{\mu_{4}}{\sigma^{4}}-3$, respectively.
Here $\mu_{n}=\langle(X-\langle X\rangle)^{n}\rangle$ and
$\sigma=\sqrt{\frac{\Sigma(X_{i}-\langle X\rangle)^{2}}{N}}$ where the $X$
represents $C_{v}^{*}$. The skewness reflects the deviation of $C_{v}^{*}$
distribution from the Gaussian distribution, and the kurtosis describes how
close the events distribute to expectancy (a standard Gaussian distribution
has zero kurtosis). At RHIC, skewness and kurtosis analysis have been
successfully applied to net-proton fluctuation to explore possible QCD
critical point [68, 69, 70, 71, 8]. Here we use this kind of high-order
moments to specific heat analysis. Fig. 8(a) and (b), respectively, showed the
skewness and kurtosis as a function of $\sqrt{s_{NN}}$, which exhibits
obviously different at lower $\sqrt{s_{NN}}$. In the String Melting version,
the skewness and kurtosis all showed a sharp enhanced structure with the
decreasing of $\sqrt{s_{NN}}$ below 20 GeV. Meanwhile in the Default version,
the skewness and kurtosis are both close to zero and present an independence
of $\sqrt{s_{NN}}$. Again, the parton interaction demonstrates its important
role for $C_{v}^{*}$ fluctuation.

Figure 8: The energy dependence of kurtosis (a) and skewness (b) of
$C_{v}^{*}$ for $\pi^{+}$ system from the String Melting AMPT (red circle) is
compared to results from the Default AMPT version (blue square) as well as
random mixed statistical fluctuations (red open circle and blue open square).

Before discussing the statistical properties obtained from $C_{v}^{*}$, we
need to take care statistical fluctuation. As Ref. [34] shows, the
contribution of statistical fluctuation in $T_{eff}$ distribution can be
written as $\left(\Delta T_{eff}\right)^{2}=\left(\Delta
T^{dyn}_{eff}\right)^{2}+\left(\Delta T^{sta}_{eff}\right)^{2}$ . In Ref.
[34], the statistical fluctuation was extracted by randomly mixing data from
experiments, and then the width of the randomly mixed temperature distribution
was removed. Here we randomly mixed each particle from a single event while
kept the multiplicity in each mixed event is close to the one in our merged
events at about 1000, then we extracted skewness and kurtosis from the mixed
events to compare with the results from two versions of AMPT.

In Figure 8(a) and (b), we can clearly find the results from the randomly
mixed events and the Default version AMPT are close to each other, which means
the energy dependence of skewness and kurtosis in the Default version without
partonic process are mainly contributed by statistical fluctuation. Meanwhile,
at lower $\sqrt{s_{NN}}$ the String Melting version showed clearly an
enhancement structure far from both the Default version and random mixed one.
This comparison demonstrates that the unique behavior of 3rd and 4th order
moment from the String Melting version at lower $\sqrt{s_{NN}}$ are mainly
from the dynamical fluctuation, which reflected that the special properties of
thermodynamic quantities indeed survived throughout partonic and hadronization
processes – even partly. Besides we can also compare the Default result and
the random mixed one, and see only slight difference. The slight difference
shows us that the self-correlation of merged events can slightly suppress
kurtosis and skewness at low $\sqrt{s_{NN}}$.

The difference between two versions, namely with or without partonic
interaction process, shows different thermal properties. When $\sqrt{s_{NN}}$
decreases, skewness of $C_{v}^{*}$ from the String Melting version increases
rapidly to a high positive value. That means, compared to the default version,
more events with higher $C_{v}^{*}$ distributed on the right side of mean
value as illustrated in Fig. 6 (a). As mentioned in this section, the real
$C_{v}^{T_{0}}$ obtained from ${}^{\prime}\xi^{\prime}$ was depressed by
averaging over the events, but the positive skewness indicated that parton
phase in the String Melting version contributes higher value of
$C_{v}^{T_{0}}$ than that in the Default AMPT version. The higher value
indicates that, in the String Melting version of AMPT, the heat capacity of
$\pi^{+}$ system behaves more similarly to prediction in other models or
theories [34, 33] than the Default version. Here the positive skewness is
explained as a signal for formation of $\pi^{+}$ system closer to the real
$C_{v}^{T_{0}}$ state, instead of explained as overall enhancement of
$C_{v}^{T_{0}}$ in $\pi^{+}$ system. Because if it is the latter case, one
should find the event-by-event distribution of $C_{v}^{*}$ will move totally
into a higher range and still follow the standard Gaussian distribution.

For the distribution of kurtosis of $C_{v}^{*}$, it is close to zero at high
$\sqrt{s_{NN}}$ and rose rapidly to large positive values below 20 GeV in the
String Melting version. Meanwhile we can see that the kurtosis from the
Default version keeps close to zero at both high and low $\sqrt{s_{NN}}$. The
low kurtosis in the Default version indicates that if we choose AMPT with
hadron gas phase instead of parton phase, the $C_{v}^{*}$ distribution shows
similarity to the standard Gaussian distribution. The sharp enhancement in the
SM version shows that the $\pi^{+}$ system formed after parton phase drives
events to a distribution with much smaller deviation ($\sigma$) than the one
formed from hadron gas. Comparing to Eq. (5) it can be found that though
$\sigma$ for $C_{v}^{*}$ and $T_{kin}$ or $T_{eff}$ are different, they both
reflect the significant change of event-by-event fluctuation of the system. As
in results from other models or data the enhancement structure of $C_{v}$ at
lower $\sqrt{s_{NN}}$ demonstrates that the fluctuation of temperature extends
to system size of all events in ensemble, and the enhancement structure of
$C_{v}^{*}$’s kurtosis indicates that the fluctuation of $C_{v}^{*}$ behaves
similarly [11, 12, 16]. The results also illustrate that the properties of
charged particles formed via parton phase are significantly different from
those experienced no partonic process. We explain the difference as a result
of partonic process, the introduction of parton provides different degrees of
freedom in hot dense matter. However, introducing parton only is not the same
as the contribution from critical phase transition, that is why we can not
extract the divergence of $C_{V}$ in the frame. But fortunately, the
significant rising trend of skewness and kurtosis in String Melting version at
lower $\sqrt{s_{NN}}$ shows us that the effect of partonic process can survive
hadronization. And further skewness and kurtosis imply us that change for the
degrees of freedom resulted from parton phase can be measured in form of high-
order moments of $C_{v}^{*}$. These observables show clear energy dependence
than $C_{V}$ from temperature fluctuation since we can obtain $C_{v}^{*}$
event-by-event distribution as many as we need. For temperature fluctuation,
however, we can only extract one $C_{V}$. In this context, one can clearly
reveal the effect from parton process by investigating $C_{v}^{*}$ and its
high-order moments. In one word, the energy dependence of $C_{v}^{*}$ together
with its skewness and kurtosis could be taken as a potential and verifiable
signal to investigate the nature around the critical point of hot dense matter
created in heavy-ion collisions.

##  V Summary

In this work, we studied the energy dependence of specific heat in a framework
of AMPT model. Specific heat was obtained via event-by-event temperature
fluctuation for $\pi^{+}$ system from the AMPT model and compared with those
from the STAR’s data [34] and other models [33], the difference among various
models and data was analyzed and the reasons were presented. Based on a few
assumptions, we derived a new quantity, namely $C_{v}^{*}$, to describe
specific heat from kinetic quantities of characteristic event. By using the
AMPT model, the $\sqrt{s_{NN}}$ dependence of $\langle C_{v}^{*}\rangle$ was
presented. This work found that the $\langle C_{v}^{*}\rangle$ behaves closer
to results from the HRG and STAR’s data than $\langle C_{v}\rangle$ from
temperature fluctuation in the AMPT frame. But at low $\sqrt{s_{NN}}$,
comparing to results from the HRG or STAR’s data, our results for $\langle
C_{v}^{*}\rangle$ still seem too small to show a clearer energy dependence.
The too small values of $\langle C_{v}^{*}\rangle$ are not the ideal results
but within the expectation.

In terms of definition in Eq. 13, our original goal of the extraction of heat
capacity should be the $\xi$ in the Mean Value Theorem for Integrals, while we
can only use the mean values of $\beta_{p_{T}}$ and $T_{eff}$ as alternative
substitutes, because the real $\xi$ is unsolvable. As a result of using the
mean value, $\langle C_{v}^{*}\rangle$ shows a visible energy dependent trend.
In other words, if one can find the event $\xi$ and use the parameters in this
event $\xi$ to calculate $C_{v}^{T_{0}}$, we would see the same energy
dependence of system as the results in experiments.

We further gave the event-by-event distribution of $C_{v}^{*}$ and then
obtained skewness and kurtosis of the $C_{v}^{*}$ distribution. To ensure the
calculation reasonable, it was additionally compared for the variance of
different terms and linear correlation coefficient. We compared parameters
like skewness and kurtosis in two different versions of AMPT.

The results demonstrate that the energy dependences of kewness from the String
Melting version show clear difference from the Default version, and similar
energy dependence can be seen in kurtosis. The above phenomenon indicates that
partonic process drives more $\pi^{+}$ events emerge in the higher specific
heat state. And from kurtosis it can be seen that the behavior of fluctuation
and deviation of $C_{v}^{*}$ in the String Melting version showed clearer
energy dependence, which is similar to the behavior of fluctuation of
temperature described by $C_{v}$ from results of Refs. [14, 36, 34, 33].

In future, we expect more experimental and theoretical investigation by using
heat capacity to search for the QCD phase transition of the hot dense matter
created in heavy-ion collisions.

###### Acknowledgements.

This work was supported in part by the National Natural Science Foundation of
China under contract Nos. 11890710, 11890714, 11875066, 11925502, 11961141003
and 12147101, the Strategic Priority Research Program of CAS under Grant No.
XDB34000000, National Key R&D Program of China under Grant No. 2018YFE0104600
and 2016YFE0100900, and by Guangdong Major Project of Basic and Applied Basic
Research No. 2020B0301030008.

## References

  * Shuryak [2017] E. Shuryak, Rev. Mod. Phys. 89, 035001 (2017).
  * Bzdak  _et al._ [2020] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanov, and N. Xu, Phys. Rep. 853, 1 (2020).
  * Braun-Munzinger  _et al._ [2016] P. Braun-Munzinger, V. Koch, T. Schäfer, and J. Stachel, Phys. Rep. 621, 76 (2016).
  * Chen  _et al._ [2018] J. Chen, D. Keane, Y.-G. Ma, A. Tang, and Z. Xu, Phys. Rep. 760, 1 (2018).
  * Tang  _et al._ [2020] Z.-B. Tang, W.-M. Zha, and Y.-F. Zhang, Nucl. Sci. Tech. 31, 81 (2020).
  * Shen and Yan [2020] C. Shen and L. Yan, Nucl. Sci. Tech. 31, 122 (2020).
  * Mrowczynski  _et al._ [2017] S. Mrowczynski, B. Schenke, and M. Strickland, Phys. Rep. 682, 1 (2017).
  * Luo and Xu [2017] X. F. Luo and N. Xu, Nucl. Sci. Tech. 28, 112 (2017).
  * Wu  _et al._ [2021a] S. Wu, C. Shen, and H. Song, Chin. Phys. Lett. 38, 081201 (2021a).
  * Waqas  _et al._ [2020a] M. Waqas, F.-H. Liu, L.-L. Li, and H. M. Alfanda, Nucl. Sci. Tech. 31, 109 (2020a).
  * Puglisi  _et al._ [2017] A. Puglisi, A. Sarracino, and A. Vulpiani, Phys. Rep. 709-710, 1 (2017).
  * Ayala  _et al._ [2019] A. Ayala, M. Hentschinski, L. Hernández, M. Loewe, and R. Zamora, Ukr. J. Phys. 64, 665 (2019).
  * Ayala  _et al._ [2018] A. Ayala, M. Hentschinski, L. A. Hernández, M. Loewe, and R. Zamora, Phys. Rev. D 98, 114002 (2018).
  * Basu  _et al._ [2015] S. Basu, R. Chatterjee, B. K. Nandi, and T. K. Nayak, arXiv e-prints , arXiv:1504.04502 (2015), arXiv:1504.04502 [nucl-ex] .
  * Shuryak [1998] E. V. Shuryak, Phys. Lett. B 423, 9 (1998).
  * Stodolsky [1995] L. Stodolsky, Phys. Rev. Lett. 75, 1044 (1995).
  * Han  _et al._ [2020] Z. Han, B. Chen, and Y. Liu, Chin. Phys. Lett. 37, 112501 (2020).
  * Pan  _et al._ [1998] J. C. Pan, S. D. Gupta, and M. Grant, Phys. Rev. Lett. 80, 1182 (1998).
  * Ma [1999] Y. G. Ma, Phys. Rev. Lett. 83, 3617 (1999).
  * D’Agostino  _et al._ [2000] M. D’Agostino, F. Gulminelli, P. Chomaz, _et al._ , Phys. Lett. B 473, 219 (2000).
  * Natowitz  _et al._ [2002] J. B. Natowitz, K. Hagel, Y. G. Ma, _et al._ , Phys. Rev. Lett. 89, 212701 (2002).
  * Gulminelli  _et al._ [2003] F. Gulminelli, P. Chomaz, A. H. Raduta, and A. R. Raduta, Phys. Rev. Lett. 91, 202701 (2003).
  * Ma  _et al._ [2005a] Y. G. Ma, J. B. Natowitz, R. Wada, _et al._ , Phys. Rev. C 71, 054606 (2005a).
  * Ma and Ma [2018] C. W. Ma and Y.-G. Ma, Prog. Part. Nucl. Phys. 99, 120 (2018).
  * Lin  _et al._ [2019] W. Lin, P. Ren, and R. Wada, Phys. Rev. C 99, 054616 (2019).
  * Borderie and Frankland [2019] B. Borderie and J. Frankland, Prog. Part. Nucl. Phys. 105, 82 (2019).
  * Liu  _et al._ [2019] H. L. Liu, Y. G. Ma, and D. Q. Fang, Phys. Rev. C 99, 054614 (2019).
  * Deng  _et al._ [2022] X. G. Deng, P. Danielewicz, Y. G. Ma, H. Lin, and Y. X. Zhang, Phys. Rev. C 105, 064613 (2022).
  * Liu  _et al._ [2022] C. Liu, X. G. Deng, and Y. G. Ma, Nucl. Sci. Tech. 33, 52 (2022).
  * Ma  _et al._ [2004] G. L. Ma, Y. G. Ma, B. H. Sa, _et al._ , High Energy Physics and Nuclear Physics 28, 398 (2004).
  * Ma  _et al._ [2005b] Y. G. Ma, G. L. Ma, X. Z. Cai, _et al._ , J. Phys. G 31, S1179 (2005b).
  * Shi  _et al._ [2007] X. H. Shi, G. L. Ma, Y. G. Ma, X. Z. Cai, and J. H. Chen, Int. J. Mod. Phys. E 16, 1912 (2007).
  * Sa  _et al._ [2007] B.-H. Sa, X.-M. Li, S.-Y. Hu, S.-P. Li, J. Feng, and D.-M. Zhou, Phys. Rev. C 75, 054912 (2007).
  * Basu  _et al._ [2016] S. Basu, S. Chatterjee, R. Chatterjee, T. K. Nayak, and B. K. Nandi, Phys. Rev. C 94, 044901 (2016).
  * Jeon and Koch [2004] S. Jeon and V. Koch, Event by event fluctuations, in  _Quark-gluon plasma. Vol. 3_, edited by R. Hwa and X. Wang (2004) pp. 430–490, arXiv:hep-ph/0304012 .
  * Deb  _et al._ [2021a] S. Deb, G. Sarwar, R. Sahoo, and J.-e. Alam, Eur. Phys. J. A 57, 195 (2021a).
  * Lin  _et al._ [2005] Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, Phys. Rev. C 72, 064901 (2005).
  * Ma and Lin [2016] G.-L. Ma and Z.-W. Lin, Phys. Rev. C 93, 054911 (2016).
  * Lin and Zheng [2021] Z.-W. Lin and L. Zheng, Nucl. Sci. Tech. 32, 113 (2021).
  * Wang and Gyulassy [1991] X.-N. Wang and M. Gyulassy, Phys. Rev. D 44, 3501 (1991).
  * Gyulassy and Wang [1994] M. Gyulassy and X.-N. Wang, Comp. Phys. Comm. 83, 307 (1994).
  * Zhang [1998] B. Zhang, Comp. Phys. Comm. 109, 193 (1998).
  * Li and Ko [1995] B.-A. Li and C. M. Ko, Phys. Rev. C 52, 2037 (1995).
  * Jin  _et al._ [2018] X.-H. Jin, J.-H. Chen, Y.-G. Ma, _et al._ , Nucl. Sci. Tech. 29, 54 (2018).
  * Wang  _et al._ [2019] H. Wang, J. H. Chen, Y. G. Ma, _et al._ , Nucl. Sci. Tech. 30, 185 (2019).
  * Lin [2014] Z.-W. Lin, Phys. Rev. C 90, 014904 (2014).
  * Mukherjee  _et al._ [2018] M. Mukherjee, S. Basu, A. Chatterjee, S. Chatterjee, S. P. Adhya, S. Thakur, and T. K. Nayak, Phys. Lett. B 784, 1 (2018).
  * Shen and Li [2020] C. Shen and Y. Li, Nucl. Sci. Tech. 31, 122 (2020).
  * Zhao  _et al._ [2019a] X.-L. Zhao, G.-L. Ma, and Y.-G. Ma, Phys. Rev. C 99, 034903 (2019a).
  * Gao  _et al._ [2020] J.-H. Gao, G.-L. Ma, S. Pu, and Q. Wang, Nucl. Sci. Tech. 31, 90 (2020).
  * Liu and Huang [2020] Y.-C. Liu and X.-G. Huang, Nucl. Sci. Tech. 31, 56 (2020).
  * Wang  _et al._ [2021] C.-Z. Wang, W.-Y. Wu, Q.-Y. Shou, G.-L. Ma, Y.-G. Ma, and S. Zhang, Phys. Lett. B 820, 136580 (2021).
  * Zhao  _et al._ [2019b] X. L. Zhao, G. L. Ma, and Y. G. Ma, Phys. Lett. B 792, 413 (2019b).
  * Wu  _et al._ [2021b] W.-Y. Wu, C.-Z. Wang, Q.-Y. Shou, Y.-G. Ma, and L. Zheng, Phys. Rev. C 103, 034906 (2021b).
  * Adamczyk  _et al._ [2017] L. Adamczyk _et al._ (STAR Collaboration), Phys. Rev. C 96, 044904 (2017).
  * Waqas  _et al._ [2020b] M. Waqas, F.-H. Liu, and Z. Wazir, Adv. High Energy Phys. 2020, 8198126 (2020b).
  * Retière and Lisa [2004] F. Retière and M. A. Lisa, Phys. Rev. C 70, 044907 (2004).
  * Chen  _et al._ [2021] J. Chen, J. Deng, Z. Tang, Z. Xu, and L. Yi, Phys. Rev. C 104, 034901 (2021), arXiv:2012.02986 [nucl-th] .
  * Su  _et al._ [2021] Y. Su, Y. J. Sun, Y. F. Zhang, and X. L. Chen, Nucl. Sci. Tech. 32, 108 (2021).
  * Deb  _et al._ [2021b] S. Deb, G. Sarwar, and R. Sahoo, The European Physical Journal A 57, 195 (2021b).
  * Landau and Lifshitz [1980] L. D. Landau and E. M. Lifshitz,  _Course of Theoretical Physics_ (1980).
  * Adams  _et al._ [2004] J. Adams _et al._ (STAR Collaboration), Phys. Rev. C 70, 054907 (2004).
  * Abelev  _et al._ [2009] B. I. Abelev _et al._ (STAR Collaboration), Phys. Rev. C 79, 034909 (2009).
  * Adam  _et al._ [2016] J. Adam _et al._ (ALICE Collaboration), Phys. Rev. C 93, 024917 (2016).
  * Adams  _et al._ [2005] J. Adams _et al._ , Nucl. Phys. A 757, 102 (2005), first Three Years of Operation of RHIC.
  * Gardim  _et al._ [2020] F. G. Gardim, G. Giacalone, M. Luzum, and J.-Y. Ollitrault, Nature Phys. 16, 615 (2020).
  * Adamczyk  _et al._ [2013] L. Adamczyk _et al._ (STAR Collaboration), Phys. Rev. C 87, 064902 (2013).
  * Stephanov [2011] M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011).
  * Aggarwal  _et al._ [2010] M. M. Aggarwal _et al._ (STAR Collaboration), Phys. Rev. Lett. 105, 022302 (2010).
  * Adamczyk  _et al._ [2014] L. Adamczyk _et al._ (STAR Collaboration), Phys. Rev. Lett. 112, 032302 (2014).
  * Adam  _et al._ [2011] J. Adam _et al._ (STAR Collaboration), Phys. Rev. Lett. 126, 092301 (2011).