Beam energy dependence of transverse momentum distribution and elliptic flow in Au-Au collisions using HYDJET++ model

Heavy Ion collision experiments at RHIC and LHC have a common objective to study an exotic state of matter called Quark Gluon Plasma(QGP) formed at extreme temperature and density, which is strong enough to de-confine the quarks and gluons [1, 2, 3]. One of the plausible approaches to study such a de-confined medium is to vary the collision energy and study its signals [4] like elliptic flow [5], jet quenching[6], strangeness abundance [7, 8], $j/\psi$ suppression [9], etc. Such an approach can be used for scanning the phase diagram of QCD as both temperature and chemical potential can be changed by varying collision energy [10, 11, 12]. Another advantage of using such a method is the possibility of finding the QCD critical point, which marks the end of the first-order phase transition [13, 14]. The critical point can be identified by enhanced fluctuations in conserved quantities like baryon number, charge, strangeness, etc [15, 16, 17].

An experimental approach to study the collision systems at wide range of beam energies was made in the RHIC Beam Energy Scan(BES) program [18, 19, 20] which explored thirteen beam energies from the top RHIC energy of $\sqrt{s_{NN}}=$ 200 GeV to the lowest energy of 3 GeV. A series of experiments under the BES program tried to find the phase boundary and critical point by looking for the signals of QGP at different collision energies. The BES program generated a large volume of experimental data for different important observables like $p_{T}$ spectra, $v_{2}$, higher moments, etc, for different collision energies. The highest energy of the BES program (i.e. 200 GeV) has been studied widely, while more phenomenological studies need to be performed for these smaller systems to study their behavior and mechanism of particle production. Furthermore, a qualitative difference between experimental data and theoretical model calculations or a deviation from the trend of model parameters at a particular beam energy could suggest a change in the particle production mechanism. Hence, the beam energy dependence of different signals and the model parameters need to be studied using a suitable theoretical model.

Some theoretical attempts were made before to predict the beam energy dependence of signals like $v_{2}$ using models like AMPT, UrQMD, etc [21, 22]. A hydrodynamical treatment can be used to get a new insight into this problem. The theoretical model we chose to employ is the HYDJET++ model. It gives the option to study the collision systems as a combination of soft QCD processes (HYDrodynamics) and hard QCD processes (JETs). Our main goal is to study the bulk properties of the QGP medium produced in lower collision energies such as $p_{T}$ spectra, particle multiplicities, and elliptic flow, which are best described by the Monte Carlo event generators. The HYDJET++ model has been successful in reproducing the experimental data at LHC energies and the top RHIC energy [23, 24, 25]. One of the main distinctions between the top RHIC energy, the LHC energies, and these lower energies is the large values of the chemical potential of these smaller systems. The model needs to be tested at lower collision energies to check the validity and implications of such treatment.

The paper is organized as follows. In section II, we have discussed the HYDJET++ framework. Section III is divided into three sub-sections. In sub-section A, we have discussed the particle ratios and freeze-out parameters. In sub-section B, we have shown the $p_{T}$ spectra at different collision energies and discussed the contribution from hard and soft processes. In sub-section C, we have shown the elliptic flow $v_{2}$ as a function of $p_{T}$ at different collision energies. Finally, in section IV, we have summarised our work.

HYDJET++ is a Monte Carlo Heavy Ion collision event generator that produces particles as a combination of hard and soft processes [24, 26]. The soft part does a detailed treatment of thermal particle production, collective flow, etc. While the hard part deals with jet production, collisional and radiative energy loss, nuclear shadowing, etc.

The hard part of the model is generated by PYQUEN. The initial patron spectra and jet production vertices from binary N-N sub-collisions are generated by PYTHIA 6 [27]. The jets produced in Au-Au collisions are a combination of jets produced from these binary sub-collisions. PYQUEN calculates the rescattering path of partons in the dense QGP medium and the radiative and collisional energy loss. The final hadronization is done according to the Lund string model. The nuclear shadowing for each impact parameter is done according to the Glauber-Gribov theory [28].

The soft part of the HYDJET++ is generated on chemical and thermal freeze-out hypersurfaces with $T_{th}$ and $T_{ch}$ as the input parameters. A detailed physics framework can be found in the corresponding papers [29, 30]. The momentum distribution of hadrons is given by the equilibrium distribution function:

Where $p^{*0}$ is the hadron energy in the fluid rest frame, $g_{i}=2J_{i}+1$ is the spin degeneracy, $\gamma^{s}$ is the strangeness suppression factor, and $\pm$ takes account of the quantum statics of a fermion or boson. The particle density at a given $T_{ch}$ and $\mu$ is given by

where $m_{i}$ is the particle mass and $K_{2}$ is the modified Bessel’s function. Chemical potential $\mu_{i}$ of a particle depends on baryon number B, strangeness S, electric charge(Isospin) Q, charm C, etc. However, for this work, we have considered only baryon number, strangeness, and i.e. $\mu_{i}=B_{i}\mu_{B}+S_{i}\mu_{s}+Q_{i}\mu_{Q}$. The baryonic chemical potential and chemical freeze-out temperature for different beam energies are chosen according to the following equation [12].

where, a=0.166$\pm$0.002 GeV, b=0.139$\pm$0.016 Ge$V^{-}1$,
c=0.053$\pm$0.021 Ge$V^{-}3$,d=1.308$\pm$0.028 GeV, e=0.273$\pm$0.008 Ge$V^{-}1$.
The values of $\mu_{s}$ are chosen to get the best match with the experimentally observed $k^{-}/k^{+}$, $k^{+}/\pi^{+}$ and $k^{-}/\pi^{-}$ ratios and $p_{T}$ spectra simultaneously. The value of isospin chemical potential is kept constant for all beam energies ($\mu_{Q}=-0.001$ GeV). The particle numbers are fixed after the chemical freeze-out. However, there are some corrections in the number of produced particles due to the decay of short-lived particles. To estimate the particle number densities at thermal freeze-out, the particle ratios at chemical freeze-out are assumed to be the same as thermal freeze-out. The effective pion chemical potential at thermal freeze-out is defined as follows,

where, $\mu_{\pi^{+}}^{eff}$ is the effective chemical potential of $\pi^{+}$ at thermal freeze-out. The values of $T_{th}$ and $\mu_{\pi^{+}}^{eff}$ are chosen to get the best match with $p_{T}$ distribution measured in RHIC BES [19]. Assuming Boltzmann approximation in eq.  (2) for particles heavier than pions [31], from eq.  (2) and eq.  (4) the chemical potential of $i^{th}$ hadron species at thermal freeze-out:

The longitudinal motion is dominant in the relativistic heavy ion collisions. Hence, the fluid flow four-velocity $\{u^{0}(x),\vec{u}(x)\}$ is parameterized in terms of longitudinal (z) and transverse ($r_{\perp}$) flow rapidities

Where $v_{\perp}$ is the magnitude of the transverse component of flow three-velocity $\vec{v}=\{v_{\perp},v_{z}\}=\{v_{\perp}cos\phi^{u},v_{\perp}sin\phi^{u},v_{z}\}$, i.e,

$u_{\perp}=\gamma v_{\perp}=cosh\eta^{u}\gamma_{\perp}v_{\perp}$, $\gamma_{\perp}=cosh\rho^{u}$. Due to the presence of non-zero azimuthal anisotropy, the azimuthal angle of the fluid velocity vector is not identical to the spatial azimuthal angle, i.e.

Where,

The parameter $\delta(b)$ is the azimuthal momentum anisotropy, a more detailed description of $\delta(b)$ is given later in the text. The transverse flow rapidity $\rho^{u}$ is related to $\rho^{{}^{\prime}u}$ as follows:

We have used a simple rapidity linear profile for $\rho^{{}^{\prime}u}$ i.e

where, r is the radial co-ordinate, $\rho_{u}^{max}(b=0)$ is the maximum transverse flow rapidity in central collisions. $R_{f}(b)$ is the mean-square radius of the hadron emission region.

In non-central collisions, the shape of the emission region in the x-y plane can be approximated by an ellipse with semi-major axis $R_{y}$ and semi-minor axis $R_{x}$ [30]. The magnitude of the $R_{f}=\sqrt{(R_{x}^{2}+R_{y}^{2})/2}$. $R_{x}(b)=R_{f}\sqrt{1-\epsilon(b)}$ and $R_{y}(b)=R_{f}\sqrt{1+\epsilon(b)}$, where $\epsilon(b)=\frac{R_{y}^{2}-R_{x}^{2}}{R_{y}^{2}+R_{x}^{2}}$. From the ellipse equation follows the dependence of the transverse radius of the hadron emission region on $\phi$

Here, $\epsilon$ is the spatial anisotropy.

The effective volume of the hadron emission region for hypersurface of proper time $\tau$ is calculated as follows:

Where $(n_{\mu}n^{\mu})=cosh\rho_{u}\sqrt{1+\delta(b)tanh^{2}\rho_{u}cos2\phi}$ and $f(\eta)$ is the longitudinal rapidity profile (assumed uniform for this work) and $\eta_{max}$ is the maximal longitudinal flow rapidity.

Elliptic flow is a result of pressure gradient in the overlap region. The expansion would be stronger along the reaction plane where the pressure gradient is larger. Due to this expansion, the initial spatial anisotropy keeps changing [32]. Since the model does not trace the evolution of spatial anisotropy, we have only used the spatial anisotropy at thermal freeze-out. The initial spatial anisotropy depends on the impact parameter and nuclear geometry. Hence, we have assumed a simple scaling between the eccentricity of the elliptical overlap region($\epsilon_{0}=b/2R_{A}$) and the spatial anisotropy at thermal freeze-out, i.e.

The parameter “$k$” is a scaling factor that indicates the change in azimuthal spatial anisotropy during evolution. It is a common practice to treat the parameters $\epsilon$ and $\delta$ as free parameters, but in such a treatment, the parameters lose all relation with the initial collision geometry. The relation can be established through the scaling factor “$k$”. The scaling factor relates the eccentricity of the overlap region $\epsilon_{0}$ with the azimuthal spatial anisotropy of the fireball at freeze-out $\epsilon$, which then determines the momentum anisotropy $\delta$. Hence, the scaling factor “$k$” determines both the spatial and momentum anisotropy parameters and the elliptic flow of final state hadrons. An empirical relation between the parameters estimated by the hydro-dynamical approach in the corresponding paper [5]:

Using the fact that $v_{2}\propto\epsilon_{0}$ and $\epsilon\propto\epsilon_{0}$, the relation between $\epsilon$ and $\delta$ can be calculated using eq.  (12),

The conversion of spatial anisotropy to momentum anisotropy, which is controlled by the parameter “C”, is assumed to be constant across all collision energies. In this work, we have fixed the value of $\mathrm{C=2}$, which was used at 200 GeV in the corresponding paper [24]. The value of the scaling factor “$k$” is chosen to get the best match with the experimental data.

Besides the above-synchronized formalism, the HYDJET++ model also has certain limitations. HYDJET++ model gives good results for the central and semi-central collision events, but it fails to give convincing results for the peripheral collisions. The model can only be used for symmetric collisions (i.e., both the colliding nuclei should be of the same element). It works reliably for relatively heavier nuclei with nuclear mass A$>$40. The model best describes the experimental data in the mid-rapidity region. Another limitation of the HYDJET++ model is that it should be used for $\sqrt{s_{NN}}>10$ GeV. For the limitations described above, we have chosen to limit our results up to 11.5 $\mathrm{GeV}$ and presented all our results in the mid-rapidity region.

We have simulated Au-Au collisions at center-of-mass energies 62.4 GeV, 39.0 GeV, 27.0 GeV, 19.6 GeV, and 11.5 GeV using the HYDJET++ event generator and compared our results with the experimental data from STAR and PHENIX detectors. The kinematic range of our results is kept the same as the experimental data. The rapidity range for particle ratios and $p_{T}$ spectra is $|y|\leq 0.1$, and for $v_{2}$ rapidity range is $|y|\leq$ 0.35 for PHENIX data and $|y|\leq$ 1.0 for STAR data. The exact number of events can be found in Table 1.

The values of freeze-out parameters such as $T_{ch}$,$\mu_{B}$, and $\mu_{s}$ used for this work are in agreement with the THERMUS model fit. A detailed analysis with the THERMUS model can be found in the corresponding papers [33, 19]. The exact values of input parameters used for this work can be found in Table 2.

We have performed a comprehensive study of bulk properties of medium like $p_{T}$ spectra and elliptic flow produced in Au-Au collisions at different beam energies. HYDJET++ model shows a good match with the experimentally observed particle ratios, which validates the Cleymans-Reidlich parameterization of $T_{ch}$ and $\mu_{B}$ for lower collision energies. HYDJET++ calculations of $p_{T}$ spectra are in agreement with experimental data for central and semi-central collisions. The lower energies of the BES have a higher chemical potential and a lower freeze-out, these conditions affect the $p_{T}$ distribution of p and $\bar{p}$. A large $\mu_{B}$ favors the thermal proton production while reducing the anti-proton production. The jet contribution also reduces due to a decrease in inelastic cross-section ($\sigma_{NN}^{inel}$) with beam energy. The interplay between these effects can be seen in the $p_{T}$ spectra of $\bar{p}$ where the jet contributions affect the low $p_{T}$ spectra. The model fails to describe the p and $\bar{p}$ spectra in peripheral collisions, this can be attributed to the centrality dependence of $\mu_{B}$, which is observed for these collision energies.

We have re-interpreted elliptic flow ($v_{2}$) based on the scaling of azimuthal spatial anisotropy. This approach allows us to calculate $v_{2}$ without having to track the whole hydrodynamic evolution of the system. The scaling between initial and final spatial anisotropy, referred to as “$k$” in this work, is used to calculate $v_{2}$. The approach can successfully describe the experimental data up to 11.5 GeV. It also preserves the essential aspects of $v_{2}$ such as mass ordering at both low $p_{T}$ ($v_{2}(\pi^{+})>v_{2}(k^{+})>v_{2}(p)$) and high $p_{T}$ ($v_{2}(p)>v_{2}(\bar{p})$). The value of “$k$” is found to be a function of beam energy and can be parameterized as a function of $\ln{\sqrt{s_{NN}}}$. The primary mesons are found to be more sensitive to $k$ than the baryons. A different value of “$k$” is required for strange and non-strange hadrons to better describe the data, which suggests that the relatively heavier strange quarks have smaller spatial and momentum anisotropy. However, this observation is empirical in nature, and more detailed studies are needed to have a better understanding of the transport properties of strange quarks.

The data will be made available on request.