# Proton number cumulants in a modified van der Waals hadron resonance gas

Kshitish Kumar Pradhan  Ronald Scaria  Dushmanta Sahu  Raghunath
Sahoo111Corresponding Author<EMAIL_ADDRESS>Department of Physics,
Indian Institute of Technology Indore, Simrol, Indore 453552, India

###### Abstract

An estimate of the proton number cumulants in the hadronic matter is presented
considering a van der Waals-type interaction between the constituent
particles. We argue that the attractive and repulsive parameters in the VDW
hadron resonance gas (VDWHRG) model change as functions of baryochemical
potential ($\mu_{B}$) and temperature ($T$). This, in turn, affects the
estimation of thermodynamic properties and, consequently, the conserved charge
fluctuations. We employ a simple parametrization to bring in the center-of-
mass energy ($\sqrt{s_{\rm NN}}$) dependence on temperature and baryochemical
potential and then estimate the proton number cumulants with the modified
approach. The modified van der Waals hadron resonance gas model (MVDWHRG)
explains the existing experimental data very well.

##  I Introduction

Exploring the quantum chromodynamics (QCD) phase diagram is of interest to the
experimental and theoretical high-energy nuclear physics community. The
existence of two phases, the hadronic phase and the quark-gluon plasma (QGP)
phase, in this diagram, separated by a crossover transition at high
temperatures ($T$) and low baryochemical potential ($\mu_{B}$), is now well
accepted based on the lattice QCD (lQCD) calculations and other
phenomenological models Bazavov:2017dus . On the other hand, it is expected
that a first-order phase transition exist between the hadronic and partonic
matter at low $T$ and high $\mu_{B}$ regime. These phase transition lines meet
at the hypothesized critical end-point (CEP). However, the applicability of
lQCD breaks down at high $\mu_{B}$ due to the fermion sign problem
Borsanyi:2013bia ; HotQCD:2014kol . Thus, to probe the CEP, which is a second-
order phase transition Stephanov:1998dy ; Stephanov:1999zu , one needs to take
the help of various phenomenological models and experiments.

It has been theorized Stephanov:1998dy ; Stephanov:1999zu ; Stephanov:2008qz ;
Hatta:2003wn  that the higher-order cumulants of conserved charges, especially
baryon number, are sensitive to the second-order phase transition. This is
because the higher-order cumulants scale with the corresponding powers of the
correlation length; thus, a limited yet finite increase of the correlation
length due to the critical slowing down can still be observed. Moreover, as
the baryon number is a conserved quantity, its fluctuations are mostly
unmodified by the final-state interactions in the hadronic phase. Some
theoretical predictions demonstrate a non-monotonic behavior of the kurtosis
of the baryon number distributions as a function of the collision energy,
given that the chemical freeze-out happens close to the CEP Stephanov:2011pb .
In an experiment, the proton number is taken as a proxy for the baryon number
given that the proton is the lightest and hence the most abundant baryon in
the hadronic matter formed after a collision. An exploration towards higher
$\mu_{B}$ by varying the collision energies in experiments to determine the
possible existence of the CEP on the phase diagram is being carried out by the
Beam Energy Scan II (BES-II) program at the Relativistic Heavy Ion Collider
(RHIC).

On the phenomenological domain, an alternative to the lQCD approach at low $T$
(up to 150 MeV) is the ideal Hadron Resonance Gas (IHRG) model. The IHRG model
has been observed to agree with the lQCD thermodynamic results for
temperatures in this range at zero $\mu_{B}$ Bellwied:2013cta ; HotQCD:2012fhj
; Bellwied:2017ttj ; Borsanyi:2011sw . It can also successfully describe the
particle ratios in high-energy collisions. Unlike lQCD, the IHRG model works
very well even at very high $\mu_{\rm B}$ regime. The IHRG model can thus be
an excellent alternative to studying the low $T$ and high $\mu_{B}$ regime.

The limiting condition of $T\leq$ 150 MeV for IHRG may be ignored while
studying the thermodynamic properties. However, this must be addressed while
estimating higher order charge fluctuations where the deviations from the lQCD
estimates are substantial Bazavov:2013dta ; Bazavov:2017dus ; Borsanyi:2018grb
; Schmidt:2012ka . Recently, much focus has been diverted towards an
interacting hadron resonance gas model as they extend the region of agreement
with lQCD data. The excluded volume hadron resonance gas (EVHRG) model assumes
an eigen volume parameter for the hadrons, which mimics a repulsive
interaction between the hadrons Andronic:2012ut ; Vovchenko:2014pka .
Meanwhile, the mean-field hadron resonance gas (MFHRG) model introduces a
repulsive interaction through an interaction potential in the hadronic medium
Kapusta:1982qd ; Olive:1980dy . There are also various other improvements to
the HRG model in literature, such as the Lorentz modified excluded volume
hadron resonance gas (LMEVHRG) model Pal:2020ink , where the hadrons are
treated as Lorentz contracted, and the effective thermal mass hadron resonance
gas (THRG) model Zhang:2019uct , where the hadrons gain effective mass with
temperature. However, the most successful improvement which explains the lQCD
results is the van der Waals hadron resonance gas (VDWHRG) model. This model
assumes a van der Waals-type interaction among the hadrons, with long-range
attractive and short-range repulsive interactions between the hadrons
Vovchenko:2015xja ; Vovchenko:2015pya ; Vovchenko:2015vxa . The VDWHRG
effectively explains the lQCD data up to $T\simeq 180$ MeV Vovchenko:2016rkn ;
Vovchenko:2017cbu . It may thus be inferred that van der Waals interaction
does play a crucial role in hadronic systems at high temperatures. In
addition, VDWHRG contains the possibility of a first-order liquid-gas phase
transition, which is a part of the QCD phase diagram. This model has been used
to estimate various thermodynamic and transport properties of the hadronic
matter Samanta:2017yhh ; Sarkar:2018mbk ; Pradhan:2022gbm . In addition, the
conserved charge fluctuations have also been estimated in the ambit of the
VDWHRG model Behera:2022nfn ; Vovchenko:2016rkn . However, very little focus
has been given to the high $\mu_{B}$ regime.

In ref. Vovchenko:2016rkn , the authors show that the VDWHRG model explains
the lQCD results better than HRG and EVHRG models. This is because repulsive
interactions significantly suppress thermodynamic functions in the crossover
region at zero $\mu_{B}$. In contrast, the attractive interactions result in
an enhancement of the thermodynamic functions. The combination of both these
interactions leads to better agreement of the model with lQCD results as
compared to the other models. In ref. Zhang:2019uct , the authors update the
VDWHRG model by considering the effective thermal mass of the hadrons. The
liquid-gas critical point has also been explored by taking the VDWHRG approach
in ref. Samanta:2017yhh . However, in all these studies, the authors have
assumed the VDW interactions to be present only among
(anti)baryons-(anti)baryons, and the meson-meson, meson-(anti)baryon
interactions have been excluded for the sake of simplicity. In ref.
Sarkar:2018mbk , a study on the criticality behavior has been conducted
including both (anti)baryons-(anti)baryons and meson-meson interactions. The
caveat is that in all of these studies Samanta:2017yhh ; Sarkar:2018mbk , the
attraction and repulsion parameters, $a$ and $b$, are taken as constants for
all $\mu_{B}$ and $T$, which seems to be an oversimplification.

In ref. Vovchenko:2016rkn , the authors have fixed the $a$ and $b$ parameters
by reproducing the saturation density $n_{\rm 0}$ = 0.16 fm-3 and binding
energy $E/A=16$ MeV of the ground state of nuclear matter. They obtain $a=329$
MeV and $b=3.42$ fm3 Vovchenko:2016rkn ; Vovchenko:2015pya ; Redlich:2016dpb
and predict a first-order liquid-gas phase transition in the nuclear matter
with a critical point at $T_{\rm c}\simeq 19.7$ MeV and $\mu_{\rm B}\simeq
908$ MeV. However, an alternate estimate of $a$ and $b$ may be obtained by
fitting the lQCD data to study the behavior of a system formed in ultra-
relativistic high-energy collisions Samanta:2017yhh . The current study argues
that the above assumptions are oversimplifications, and the attractive and
repulsive parameters may depend on chemical potential and temperature. Using a
chi-square minimization technique, we fit the $\mu_{\rm B}/T$ dependent lQCD
data to obtain the $a$ and $b$ parameters at different $\mu_{B}/T$. We further
obtain a functional form for $a$ and $b$ as functions of $\mu_{\rm B}/T$,
which can be used for any reasonably extrapolated $\mu_{B}$ and $T$ values.
This essentially modifies the results of the VDWHRG model at a high baryon-
density region. Taking this modified VDWHRG (MVDWHRG) approach, we try to
explain the proton number fluctuation data from the STAR collaboration.

The current work is organized as follows. In section II, we briefly describe
the van der Waals hadron resonance gas model. We then discuss the modification
of the $a$ and $b$ parameters as functions of baryochemical potential. In
section III, we briefly discuss the obtained results and consequences.
Finally, we summarise our results in section IV.

##  II Formulation

The ideal HRG formalism considers hadrons to be point particles with no
interactions between them. Under this formalism, the partition function of
$i^{th}$ particle species having mass $m_{i}$ is given in a grand canonical
ensemble (GCE) as Andronic:2012ut ,

$lnZ^{id}_{i}=\pm\frac{Vg_{i}}{2\pi^{2}}\int_{0}^{\infty}p^{2}dp\
ln\\{1\pm\exp[-(E_{i}-\mu_{i})/T]\\},$ (1)

where the notations of $g_{i}$, $E_{i}=\sqrt{p^{2}+m_{i}^{2}}$ and $\mu_{i}$
are the degeneracy, energy and chemical potential of the $i^{th}$ hadron,
respectively. $\mu_{i}$ is further expanded in terms of the baryonic,
strangeness, and charge chemical potentials ($\mu_{B}$, $\mu_{S}$ and
$\mu_{Q}$, respectively) and the corresponding conserved numbers ($B_{i}$,
$S_{i}$ and $Q_{i}$) as,

$\mu_{i}=B_{i}\mu_{B}+S_{i}\mu_{S}+Q_{i}\mu_{Q},$ (2)

Pressure $P_{i}$, energy density $\varepsilon_{i}$, number density $n_{i}$,
and entropy density $s_{i}$ in the ideal HRG formalism can now be obtained as,

$P^{id}_{i}(T,\mu_{i})=\pm\frac{Tg_{i}}{2\pi^{2}}\int_{0}^{\infty}p^{2}dp\
ln\\{1\pm\exp[-(E_{i}-\mu_{i})/T]\\}$ (3)
$\varepsilon^{id}_{i}(T,\mu_{i})=\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}\frac{E_{i}\
p^{2}dp}{\exp[(E_{i}-\mu_{i})/T]\pm 1}$ (4)
$n^{id}_{i}(T,\mu_{i})=\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}\frac{p^{2}dp}{\exp[(E_{i}-\mu_{i})/T]\pm
1}$ (5) $\displaystyle s^{id}_{i}(T,\mu_{i})=$
$\displaystyle\pm\frac{g_{i}}{2\pi^{2}}\int_{0}^{\infty}p^{2}dp\Big{[}\ln\\{1\pm\exp[-(E_{i}-\mu_{i})/T]\\}$
$\displaystyle\pm\frac{(E_{i}-\mu_{i})/T}{\exp[(E_{i}-\mu_{i})/T]\pm
1}\Big{]}.$ (6)

8

Figure 1: (Color Online) Simultaneous $\chi^{2}$ fit of scaled energy density
and pressure at various $\mu_{B}/T$ as functions of temperature. The coloured
bands are for lQCD data Bazavov:2017dus , and the solid lines are the fits to
the data.

Figure 2: (Color Online) Variation of extracted $a$ and $b$ parameters from
the $\chi^{2}$ fit to the lQCD data as functions of $\mu_{B}/T$. The solid
black lines denote the negative exponential fits to the $a$ and $b$
parameters.

The ideal HRG model can be modified to include van der Waals interactions
between particles by the introduction of the attractive and repulsive
parameters $a$ and $b$, respectively. This modifies the pressure and number
density obtained in ideal HRG iteratively as follows Samanta:2017yhh ;
Vovchenko:2015vxa ; Vovchenko:2015pya ;

$P(T,\mu)=P^{id}(T,\mu^{*})-an^{2}(T,\mu),$ (7)

where the $n(T,\mu$) is the VDW particle number density given by,

$n(T,\mu)=\frac{\sum_{i}n_{i}^{id}(T,\mu^{*})}{1+b\sum_{i}n_{i}^{id}(T,\mu^{*})}.$
(8)

Here, $i$ runs over all hadrons and $\mu^{*}$ is the modified chemical
potential given by,

$\mu^{*}=\mu-bP(T,\mu)-abn^{2}(T,\mu)+2an(T,\mu).$ (9)

It is to be noted that the repulsive parameter is usually attributed to be
related to the hardcore radius ($r$) of the particle through the relation
$b=16\pi r^{3}/3$, while no such definition is known for the attractive
parameter.

Entropy density $s(T,\mu)$ and energy density $\epsilon(T,\mu)$ in VDWHRG can
now be obtained as,

$s(T,\mu)=\frac{s^{id}(T,\mu^{*})}{1+bn^{id}(T,\mu^{*})}$ (10)
$\epsilon(T,\mu)=\frac{\sum_{i}\epsilon_{i}^{id}(T,\mu^{*})}{1+b\sum_{i}n_{i}^{id}(T,\mu^{*})}-an^{2}(T,\mu).$
(11)

Knowing the attractive and repulsive parameters, the total pressure in the
VDWHRG model can be written as Samanta:2017yhh ; Vovchenko:2015vxa ;
Vovchenko:2015pya ; Vovchenko:2017cbu ; Vovchenko:2016rkn ; Sarkar:2018mbk ,

$P(T,\mu)=P_{M}(T,\mu)+P_{B}(T,\mu)+P_{\bar{B}}(T,\mu)$ (12)

where $M$, $B$ and $\bar{B}$ represents mesons, baryons and anti-baryons,
respectively. Their contributions to pressure is now defined by,

$P_{M}(T,\mu)=\sum_{i\in M}P_{i}^{id}(T,\mu^{*M})$ (13)
$P_{B}(T,\mu)=\sum_{i\in B}P_{i}^{id}(T,\mu^{*B})-an^{2}_{B}(T,\mu)$ (14)
$P_{\bar{B}}(T,\mu)=\sum_{i\in\bar{B}}P_{i}^{id}(T,\mu^{*\bar{B}})-an^{2}_{\bar{B}}(T,\mu).$
(15)

The excluded volume correction modifies the meson chemical potential as
$\mu^{*M}$, while VDW interactions give modified chemical potentials of
baryons and anti-baryons as $\mu^{*B}$ and $\mu^{*\bar{B}}$ Sarkar:2018mbk .
For the simple case of vanishing electric charge and strangeness chemical
potentials Bazavov:2017dus , $\mu_{Q}=\mu_{S}=0$, these modified chemical
potentials are obtained from Eq. 2 and Eq. 9 as;

$\mu^{*M}=-bP_{M}(T,\mu)$ (16)
$\mu^{*B(\bar{B})}=\mu_{B(\bar{B})}-bP_{B(\bar{B})}(T,\mu)-abn^{2}_{B(\bar{B})}+2an_{B(\bar{B})}$
(17)

where $n_{M}$, $n_{B}$ and $n_{\bar{B}}$ denotes modified number densities of
each type of particle given by

$n_{M}(T,\mu)=\frac{\sum_{i\in M}n_{i}^{id}(T,\mu^{*M})}{1+b\sum_{i\in
M}n_{i}^{id}(T,\mu^{*M})}$ (18) $n_{B(\bar{B})}(T,\mu)=\frac{\sum_{i\in
B(\bar{B})}n_{i}^{id}(T,\mu^{*B(\bar{B})})}{1+b\sum_{i\in
B(\bar{B})}n_{i}^{id}(T,\mu^{*B(\bar{B})})}.$ (19)

Here the summation runs over mesons ($M$), baryons ($B$), and anti-baryons
($\bar{B}$) respectively.

As already discussed, the VDW parameters $a$ and $b$ may depend upon
temperature and chemical potential. In ref. Dutra:2020qsn , the authors have
considered a density-dependent VDW model where the VDW parameters are no
longer constants but are defined as functions of density. Taking forward the
assumption of a fixed hard core radius for mesons, we have estimated the VDW
parameters for (anti)baryon-(anti)baryon interactions for different
$\mu_{B}/T$ values by fitting the corresponding lQCD energy density and
pressure Bazavov:2017dus . In order to extract the van der Waals parameters
$a$ and $b$, we use the $\chi^{2}$ minimization technique to obtain the best
fit to the available lQCD data Bazavov:2017dus  using the relation
Samanta:2017yhh ; Vovchenko:2014pka ,

$\chi^{2}=\frac{1}{N_{df}}\sum_{i,j}\frac{\big{(}X^{lQCD}_{i,j}(T_{j})-X^{model}_{i,j}(T_{j})\big{)}^{2}}{\big{(}\Delta^{lQCD}_{i,j}(T_{j})\big{)}^{2}},$
(20)

where the $X^{lQCD}_{i,j}(T_{j})$ and $X^{model}_{i,j}(T_{j})$ are the values
of the $i^{th}$ thermodynamic observable from lQCD calculations and model
calculations respectively at a given $j^{th}$ temperature.
$\Delta^{lQCD}_{i,j}$ is the error involved in the lQCD calculation. $N_{df}$
is the number of lQCD data points minus the number of fitting parameters. Fig.
1 shows the goodness of the fit (solid line) to the lQCD data (color band) for
different values of $\mu_{B}/T$ for scaled pressure and energy density. The
obtained parameters and the corresponding $\chi^{2}$ values for each case of
$\mu_{B}/T$ are listed in Table 1. The attractive and repulsive parameters,
$a$ and $b$, are now shown as a function of $\mu_{B}/T$ in fig. 2. It is
observed that both these parameters decrease as functions of $\mu_{B}/T$.
Therefore an increase in the baryon density or a decrease in the temperature
leads to the weakening of the van der Waals interaction. This behavior is
similar to what was observed in ref Dutra:2020qsn . To parameterize these two
parameters as functions of $\mu_{B}$ and $T$, the obtained $a$ and $b$ values
are fitted with a negative exponential function. Therefore the van der Waals
parameters now vary as functions of $\mu_{B}/T$, quantified by the relation,

$\begin{split}a=c_{1}\exp(c_{2}\frac{\mu_{B}}{T}),\\\
b=d_{1}\exp(d_{2}\frac{\mu_{B}}{T}).\end{split}$ (21)

where the constants $c_{1}$, $c_{2}$, $d_{1}$, and $d_{2}$ are given by,
$c_{1}$ = 1.66 $\pm$ 0.05 GeV fm3, $c_{2}$ = -0.88 $\pm$ 0.04, $d_{1}$ =
541.93 $\pm$ 15.98 GeV-3, and $d_{2}$ = -0.61 $\pm$ 0.03. As it is evident
from fig. 2, both the VDW parameters are decreasing exponentially with
$\mu_{B}/T$, therefore at high baryon density and low temperature the
interacting hadron resonance gas approaches towards the ideal gas condition.
In view of the van der Waals parameters getting modified with $\mu_{B}/T$, the
model is henceforth referred to as the modified van der Waals HRG (MVDWHRG)
model.

$\mu_{B}/T$ |  $a$ (GeV fm3) |  $r_{B}$ (fm) | $\chi^{2}$
---|---|---|---
0.0 | 1.650 $\pm$ 0.05 | 0.635 $\pm$ 0.05 | 1.06/20
1.0 | 0.786 $\pm$ 0.064 | 0.515 $\pm$ 0.05 | 0.91/20
2.0 | 0.275 $\pm$ 0.025 | 0.425 $\pm$ 0.05 | 1.88/20
2.5 | 0.150 $\pm$ 0.05 | 0.385 $\pm$ 0.15 | 3.5/20
Table 1: van der Waals parameters obtained using $\chi^{2}$ minimization
technique for different $\mu_{B}/T$ lQCD data. Here $r_{\rm B}$ is the baryon
radius.

Figure 3: (Color Online) Dimensionless ratios of net proton cumulants as
functions of center-of-mass energies. The blue, green, and red solid lines
represent results obtained with IHRG, VDWHRG and MVDWHRG respectively. The
obtained results are compared with the STAR data STAR:2021iop  represented by
black markers.

The fluctuations of conserved charges can indicate the existence of a second-
order phase transition and hence the position of the critical point in the QCD
phase diagram. Such a point is a thermodynamic singularity where the
susceptibilities diverge, and the order parameter fluctuates Stephanov:1999zu
. Therefore the susceptibilities and their ratios show non-monotonic behavior
as a function of collision energy. Fluctuations of a conserved quantity and
its higher moments can be obtained from the derivatives of the basic
thermodynamic quantity, pressure, with respect to the corresponding chemical
potential. These derivatives are known as susceptibilities or cumulants. The
net proton number fluctuations can be studied as a proxy for the net baryon
number fluctuations. The susceptibilities are given by,

$\kappa_{n}=\frac{\partial^{n}}{\partial(\mu_{B}/T)^{n}}\frac{P}{T^{4}}.$ (22)

where the $P$ is the pressure defined in Eq. 3 and with the inclusion of van
der Waals interaction, it is given by Eq. 12. From Eq. 22, the first
derivative gives the mean of the net proton number as, $M=VT^{3}\kappa_{1}$,
where $\kappa_{1}$ is simply the net proton number density given by
$(n_{p}-n_{\bar{p}})/T^{3}$. Similarly, one can define the variance (Gaussian
width), $\sigma^{2}$ along with the non-Gaussian fluctuations such as the
skewness, $S$, and the kurtosis, $\kappa$ Stephanov:2008qz ; Karsch:2010ck ;
Fukushima:2014lfa ; Vovchenko:2021kxx , so as to define the cumulant ratios
as,

$S\sigma=\frac{\kappa_{3}}{\kappa_{2}},\hskip
14.22636pt\frac{S\sigma^{3}}{M}=\frac{\kappa_{3}}{\kappa_{1}},\hskip
14.22636pt\kappa\sigma^{2}=\frac{\kappa_{4}}{\kappa_{2}}.$ (23)

We can now estimate these ratios using the MVDWHRG model. However, in order to
compare our results with the experimentally available data, we need to relate
the parameters of the hadron resonance gas, $\mu_{B}$ and $T$ with the
collision energy, $\sqrt{s_{\rm NN}}$. Such a parametrization is available in
literature Cleymans:2005xv  and is given by,

$T(\mu_{B})=q_{1}-q_{2}\mu_{B}^{2}-q_{3}\mu_{B}^{4},$ (24)
$\mu_{B}(\sqrt{s_{NN}})=\frac{q_{4}}{1+q_{5}\sqrt{s_{NN}}},$ (25)

where the fitting parameters are given as $q_{1}$ = 0.166 GeV, $q_{2}$ = 0.139
GeV-1, $q_{3}$ = 0.053 GeV-3, $q_{4}$ = 1.308 GeV, and $q_{5}$ = 0.273 GeV-1.
Note that these parameters are obtained by using a freeze-out criterion based
on the IHRG model and are not quite model dependent Cleymans:2005xv . There is
a very minimal change in the parameters when including interactions like
excluded volume (EV) interactions Tiwari:2011km  or the conventional VDW
interaction Behera:2022nfn . Therefore, we use the IHRG parameters throughout
this work.

The expression for the proton number fluctuation given in Eq. (22) is modified
from ideal HRG owing to the van der Waals interaction. The baryon chemical
potential $\mu_{B}$ is reduced to modified chemical potential $\mu_{B}^{*}$ as
defined in Eq. 15. Similarly, the pressure and number density are also
modified as in Eq. 10, and Eq. 17. This leads to a large deviation in the
higher-order cumulants in MVDWHRG from ideal HRG, as will be shown in the next
section. We have used Eq. 24 and Eq. 25 for the $\sqrt{s_{NN}}$ dependencies
of the desired ratios in Eq. 23.

##  III Results and Discussion

It is expected that the medium formed in relativistic heavy ion collisions
undergo particlization at the chemical freeze-out transforming into a hadron
resonance gas before further evolution towards the detectors. In mapping the
QCD phase diagram, conserved charge fluctuations have been proposed and
explored as a method to pinpoint the existence and location of the QCD
critical point. Net-proton fluctuations, in particular, have been seen as a
proxy for the baryon number fluctuations as they constitute the dominant part
of the baryons that can be detected. Multiple studies have tried to relate the
proton number fluctuations to the baryon number fluctuations within the ambit
of various models with methods available in the literature to relate the two
Fukushima:2014lfa ; Vovchenko:2021kxx ; Kitazawa:2011wh ; Kitazawa:2012at .

Figure 4: (Color Online) Ratios of proton number cumulants as functions of the
center-of-mass energy. The red lines and marker correspond to the proton and
the blue lines and marker correspond to the anti-proton. The results obtained
in MVDWHRG are shown by the solid line, and that of VDWHRG are by the dashed
line. The dotted lines are for the results from hydrodynamical calculations
Vovchenko:2021kxx , whereas, the solid markers are for the STAR data
STAR:2021iop .

Figure 5: (Color Online) Correlation functions as functions of center-of-mass
energy. The red lines and marker correspond to the proton and the blue lines
and marker correspond to the anti-proton. The results obtained in MVDWHRG are
shown by the solid line, and that of VDWHRG are by the dashed line. The dotted
lines are for the results from hydrodynamical calculations Vovchenko:2021kxx ,
whereas, the solid markers are for the STAR data STAR:2021iop .

In principle, the interactions among the hadrons quantify the various dynamics
of the hadronic medium. Thus, by considering the attraction and repulsion
between the hadrons, one can understand the dynamical behavior of the matter
formed in high-energy collisions. In this regard, we have incorporated the
baryochemical and temperature dependence of the van der Waals interaction
parameters in section II, thus appropriately modifying the equation of state.
Now, with the modified $a$ and $b$ parameters, we proceed to estimate the net-
proton fluctuations as a function of the center-of-mass energy.

In the Boltzmann limit, i.e. when the number densities of both $p$ and
$\bar{p}$ are both dilute, one may obtain the Skellam limit for the ratios
defined in Eq. 23,

$S\sigma={\rm tanh}~{}(\mu_{B}/T),\hskip
14.22636pt\frac{S\sigma^{3}}{M}=1,\hskip 14.22636pt\kappa\sigma^{2}=1.$ (26)

Any deviation from this limit may indicate the onset of new physics. The
presence of QCD critical point has been predicted to strongly affect these
quantities Stephanov:2008qz ; Stephanov:2011pb ; Mroczek:2020rpm . Fig.3 shows
the variation of these parameters as functions of the center-of-mass energy.
The solid black markers show the net-proton fluctuation values obtained from
the RHIC BES-I experiment STAR:2021iop . The black dashed line represents the
Skellam predictions, and the solid blue line shows ideal HRG results. The
VDWHRG values in solid green are obtained by using the constant VDW
parameters, $a$ = 329 MeV fm3 and $b$ = 3.42 fm3 Vovchenko:2015pya . All these
results are compared with the modified VDWHRG values (solid red), which are
obtained by using the functional form of the $a$ and $b$ parameters obtained
in section II. It was shown in ref. Fukushima:2014lfa  that the HRG model
deviates from the Skellam condition at high baryon densities due to the
inclusion of quantum correlations. But this alone is insufficient to explain
the BES-I data for the proton number cumulants, suggesting that these quantum
correlations are weak Kitazawa:2011wh ; Kitazawa:2012at . The authors in ref.
Vovchenko:2021kxx  used a hybrid model where further interactions were
considered, like the excluded volume effect and global baryon conservation,
which were seen to have a significant effect on these ratios. It was suggested
that a simple excluded volume component explains the higher center-of-mass
energy region well while the deviation below $\sqrt{s_{NN}}\leq 20$ GeV
signifies the importance of attractive interactions. The VDWHRG model provides
the opportunity to test this but fails in explaining the data, while the
MVDWHRG model described in this paper shows a good agreement with the obtained
net-proton number fluctuations from BES-I. It can be observed that the
deviation from the VDWHRG and MVDWHRG is lesser near the high center-of-mass
energies (low $\mu_{B}$). In contrast, the deviation increases with decreasing
the center-of-mass energies (high $\mu_{B}$). The interesting thing to note
here is that at low center-of-mass energies (high $\mu_{B}$), the MVDWHRG
matches with the ideal HRG. This is because at low center-of-mass energy,
which corresponds to high values of baryon density, $\mu_{B}$, both of the VDW
parameters weaken, approaching zero. Hence the system behaves more like an
ideal one.

Further, one may explore the medium by studying the cumulants of protons and
anti-protons separately. Correlation functions, $C_{n}$, may be constructed as
a linear combination of ordinary cumulants Vovchenko:2021kxx  as given below,

$C_{1}=\kappa_{1}$ (27) $C_{2}=-\kappa_{1}+\kappa_{2}$ (28)
$C_{3}=2\kappa_{1}-3\kappa_{2}+\kappa_{3}$ (29)
$C_{4}=-6\kappa_{1}+11\kappa_{2}-6\kappa_{3}+\kappa_{1}.$ (30)

These values approach zero for $n>1$ in the Poissonian limit, indicating
uncorrelated particle production. Both $\kappa_{n}$ and $C_{n}$ were explored
at STAR STAR:2021iop  as functions of the center-of-mass energy. The
comparison of our calculations of these quantities compared with the
experimental results is shown in Fig. 4 and 5. The red and cyan colors
represent the results for proton and anti-protons, respectively. The
deviations observed from zero indicate the existence of physics beyond the
standard uncorrelated hadron gas scenario. The normalized quantities,
$\kappa_{n}/\kappa_{1}-1$ and $C_{n}/C_{1}$, are well explained by the MVDWHRG
model as compared to the standard VDWHRG model as can be seen from the solid
and dashed curves. The results from viscous hydrodynamic simulations that
include the combined effect of the excluded volume and baryon conservation
were shown to have significant contributions at $\sqrt{s_{NN}}\geq 20$ GeV
region as observed in ref. Vovchenko:2021kxx . The change of interaction from
repulsive to attractive was suggested as a solution to the deviation observed
at lower energies. The inclusion of attractive interactions, which competes
with the repulsive interaction through the functional form of $a$ and $b$ is
observed to explain the effects being observed in data. These results suggest
that using constant $a$ and $b$ obtained by fitting lQCD results at
$\mu_{B}/T=0$ is not sufficient to describe the system at a non-zero $\mu_{B}$
adequately. It may also be noted that only a qualitative description of the
anti-proton trend is obtained. But, even then, our results are better in
agreement with the data as compared with the results obtained in
hydrodynamical models Vovchenko:2021kxx  or the URQMD model Xu:2016qjd ;
He:2017zpg  which was used for comparison with data in the STAR experiment
STAR:2021iop .

Finally, In Fig. 6, we present our predictions for the net proton
hyperkurtosis, $\kappa_{6}/\kappa_{2}$. The result from the ideal HRG shown by
the solid blue line is decreasing towards the low center-of-mass energy but
remains positive throughout the selected values of $\sqrt{s_{NN}}$. However,
the red solid line that represents the results from the MVDWHRG remains
negative at the high center-of-mass energy and becomes positive towards
$\sqrt{s_{NN}}\leq 5$ and matches with ideal HRG. The prediction for
$\kappa_{6}/\kappa_{2}$ from ref. Vovchenko:2021kxx  is shown for comparison.

Figure 6: (Color Online) Hyperkurtosis as a function of center-of-mass energy.
The ideal HRG is shown by a blue solid curve and the results obtained from
MVDWHRG is represented by a red solid curve. The dotted-dashed line is for the
results obtained using hydrodynamic calculation Vovchenko:2021kxx .

##  IV Summary

In this paper, we present an estimate of the van der Waals (VDW) attractive
($a$) and repulsive ($b$) parameters obtained from reduced $\chi^{2}$ fits to
the lQCD thermodynamic quantities at fixed $\mu_{B}/T$. The obtained values
are then parameterized as exponential functions of $\mu_{B}/T$. We have
further estimated the proton number fluctuations within this modified van der
Waals HRG (MVDWHRG) model as a function of center-of-mass energy. It is
observed that the MVDWHRG model is able to explain the proton number
fluctuation data from the RHIC BES-I experiment. We also give predictions for
higher-order cumulants and hyperkurtosis, which is yet to be measured in data.

It is to be noted that the generally expected liquid-gas phase transition, due
to VDW interactions, is absent in our model. This is due to the exponential
decrease of $a$ and $b$ as a function of $\mu_{B}/T$, thus approaching the
ideal HRG model at high $\mu_{B}/T$. But as seen from fits to the $\chi$EFT
data at $T=0$ GeV Fujimoto:2021dvn , the VDW $a$ and $b$ parameters lie above
the lowest values obtained by us on fitting lQCD results. This indicates that
the VDW parameters must reach a minimum and then rise towards higher
$\mu_{B}/T$ thus deviating from the exponential fit used in this work and
pointing towards possible critical behaviour. But due to the paucity of lQCD
results and also of $\chi$EFT data a study on the VDW parameters at these
values of $\mu_{B}/T$ is at present not possible. It would thus be worthwhile
and exciting to look for novel approaches to reconcile these first principle
models, thereby allowing a complete study of the $T-\mu_{B}$ plane.

## Acknowledgement

KP and Ronald Scaria acknowledge the financial aid from UGC and CSIR,
Government of India, respectively. The authors gratefully acknowledge the DAE-
DST, Govt. of India funding under the mega-science project – “Indian
participation in the ALICE experiment at CERN” bearing Project No.
SR/MF/PS-02/2021-IITI (E-37123).

## References

  * (1) A. Bazavov, H. T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa, S. Mukherjee, H. Ohno and P. Petreczky, et al. Phys. Rev. D 95, 054504 (2017).
  * (2) S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabo, Phys. Lett. B 730, 99 (2014).
  * (3) A. Bazavov et al. [HotQCD], Phys. Rev. D 90, 094503 (2014).
  * (4) M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998).
  * (5) M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. D 60, 114028 (1999).
  * (6) Y. Hatta and M. A. Stephanov, Phys. Rev. Lett. 91, 102003 (2003) [erratum: Phys. Rev. Lett. 91, 129901 (2003)].
  * (7) M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009).
  * (8) M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011).
  * (9) R. Bellwied, S. Borsanyi, Z. Fodor, S. D. Katz and C. Ratti, Phys. Rev. Lett. 111, 202302 (2013).
  * (10) A. Bazavov et al. [HotQCD], Phys. Rev. D 86, 034509 (2012).
  * (11) R. Bellwied, S. Borsanyi, Z. Fodor, J. Gunther, K. H. Kampert, S. D. Katz, T. Kawanai, T. G. Kovacs, S. W. Mages and A. Pasztor, et al. Nucl. Phys. A 967, 732 (2017).
  * (12) S. Borsanyi, Z. Fodor, S. D. Katz, S. Krieg, C. Ratti and K. Szabo, JHEP 01, 138 (2012).
  * (13) A. Bazavov, H. T. Ding, P. Hegde, O. Kaczmarek, F. Karsch, E. Laermann, Y. Maezawa, S. Mukherjee, H. Ohno and P. Petreczky, et al. Phys. Rev. Lett. 111, 082301 (2013).
  * (14) C. Schmidt [BNL-Bielefeld], Nucl. Phys. A 904, 865c (2013).
  * (15) S. Borsanyi, Z. Fodor, J. N. Guenther, S. K. Katz, K. K. Szabo, A. Pasztor, I. Portillo and C. Ratti, JHEP 10, 205 (2018).
  * (16) A. Andronic, P. Braun-Munzinger, J. Stachel and M. Winn, Phys. Lett. B 718, 80 (2012).
  * (17) V. Vovchenko, D. V. Anchishkin and M. I. Gorenstein, Phys. Rev. C 91, 024905 (2015).
  * (18) J. I. Kapusta and K. A. Olive, Nucl. Phys. A 408, 478 (1983).
  * (19) K. A. Olive, Nucl. Phys. B 190, 483 (1981).
  * (20) S. Pal, A. Bhattacharyya and R. Ray, Nucl. Phys. A 1010, 122177 (2021).
  * (21) H. X. Zhang, J. W. Kang and B. W. Zhang, Phys. Rev. D 101, 114033 (2020).
  * (22) V. Vovchenko, D. V. Anchishkin and M. I. Gorenstein, J. Phys. A 48, 305001 (2015).
  * (23) V. Vovchenko, D. V. Anchishkin, M. I. Gorenstein and R. V. Poberezhnyuk, Phys. Rev. C 92, 054901 (2015).
  * (24) V. Vovchenko, D. V. Anchishkin and M. I. Gorenstein, Phys. Rev. C 91, 064314 (2015).
  * (25) V. Vovchenko, M. I. Gorenstein and H. Stoecker, Phys. Rev. Lett. 118, 182301 (2017).
  * (26) V. Vovchenko, Phys. Rev. C 96, 015206 (2017).
  * (27) S. Samanta and B. Mohanty, Phys. Rev. C 97, 015201 (2018).
  * (28) N. Sarkar and P. Ghosh, Phys. Rev. C 98, 014907 (2018).
  * (29) K. K. Pradhan, D. Sahu, R. Scaria and R. Sahoo, Phys. Rev. C 107, 014910 (2023).
  * (30) S. P. Behera and D. K. Mishra, Nucl. Phys. A 1024, 122475 (2022).
  * (31) K. Redlich and K. Zalewski, Acta Phys. Polon. B 47, 1943 (2016).
  * (32) M. Dutra, B. M. Santos and O. Lourenço, J. Phys. G 47, 035101 (2020).
  * (33) F. Karsch and K. Redlich, Phys. Lett. B 695, 136 (2011).
  * (34) K. Fukushima, Phys. Rev. C 91, 044910 (2015).
  * (35) V. Vovchenko, V. Koch and C. Shen, Phys. Rev. C 105, 014904 (2022).
  * (36) J. Cleymans, H. Oeschler, K. Redlich and S. Wheaton, Phys. Rev. C 73, 034905 (2006).
  * (37) S. K. Tiwari, P. K. Srivastava and C. P. Singh, Phys. Rev. C 85, 014908 (2012).
  * (38) M. Kitazawa and M. Asakawa, Phys. Rev. C 85, 021901 (2012).
  * (39) M. Kitazawa and M. Asakawa, Phys. Rev. C 86, 024904 (2012) [E: Phys. Rev. C 86, 069902 (2012)].
  * (40) D. Mroczek, A. R. Nava Acuna, J. Noronha-Hostler, P. Parotto, C. Ratti and M. A. Stephanov, Phys. Rev. C 103, 034901 (2021).
  * (41) M. Abdallah et al. [STAR], Phys. Rev. C 104, 024902 (2021).
  * (42) J. Xu, S. Yu, F. Liu and X. Luo, Phys. Rev. C 94, 024901 (2016).
  * (43) S. He and X. Luo, Phys. Lett. B 774, 623 (2017).
  * (44) Y. Fujimoto, K. Fukushima, Y. Hidaka, A. Hiraguchi and K. Iida, Phys. Lett. B 835, 137524 (2022).