(3+1)-D viscous hydrodynamics CLVisc at finite net baryon density: identified particle spectra, anisotropic flows and flow fluctuations across BES energies

In this work, we use the Monte-Carlo Glauber model to provide initial conditions for the hydrodynamic evolution. The local entropy density $s(x,y,\eta_{s})$ and local baryon density $n(x,y,\eta_{s})$ at the initial proper time $\tau_{0}$ are obtained according to the following expressions (Denicol:2018wdp, ):

	$\displaystyle s(x,y,\eta_{s})=\frac{K}{\tau_{0}}(H^{s}_{P}(\eta_{s})s_{P}(x,y)+H^{s}_{T}(\eta_{s})s_{T}(x,y))\,,$		(1)
	$\displaystyle n(x,y,\eta_{s})=\frac{1}{\tau_{0}}(H^{n}_{P}(\eta_{s})s_{P}(x,y)+H^{n}_{T}(\eta_{s})s_{T}(x,y))\,,$		(2)

where $s_{P}(x,y)$ and $s_{T}(x,y)$ are the entropy densities in the transverse plane produced from the projectile and target nuclei, $H_{P}^{s}(\eta_{s})$, $H_{T}^{s}(\eta_{s})$, $H_{P}^{n}(\eta_{s})$ and $H_{T}^{n}(\eta_{s})$ are longitudinal envelop functions which describe the rapidity dependences for the initial entropy density and baryon number density associated with projectile and target nuclei, $K$ is the scale factor to control the magnitude of initial entropy, and $\tau_{0}$ is the initial proper time for the hydrodynamic evolution. Since the Lorentz contraction is smaller at low collision energies, the time for two nuclei to pass each other is larger.

One may estimate the overlap time as $\tau_{\rm overlap}=\frac{2R}{\sinh(y_{\rm beam})}$ , where $R$ denotes the radius of the nucleus and $y_{\rm beam}$ represents the beam rapidity. In this work, we do not include the pre-equilibrium stage, which becomes more important in lower collision energies (Dore:2020fiq, ). To compensate some effects from the pre-equilibrium evolution, we choose the initial proper time $\tau_{0}$ of hydrodynamics evolution to be larger than the overlap time $\tau_{\rm overlap}$ of two nuclei so that the system has more time to approach the local equilibrium before hydrodynamics starts.

As for the entropy density $s_{P/T}(x,y)$ in the transverse plane, it is the sum of the Gaussian smearing functions over all the participating projectile and target nucleons in the Monte-Carlo Glauber model:

$$s_{P/T}(x,y)=\sum_{i}^{P/T}\frac{1}{2\pi\sigma_{r}^{2}}\exp\frac{-(x_{i}-x)^{2}-(y_{i}-y)^{2}}{2\sigma_{r}^{2}}\,,$$

(3)

where $(x_{i},y_{i})$ is the transverse position of the participant nucleon and $\sigma_{r}$ is the transverse Gaussian smearing width (we take $\sigma_{r}=0.5$ fm). The longitudinal envelop functions $H^{s}_{P/T}$ and $H^{n}_{P/T}$ are parametrized as follows:

	$\displaystyle H^{s}_{P/T}(\eta_{s})$	$\displaystyle=\theta(\eta_{\text{max}}-|\eta_{s}|)\left(1\pm\frac{\eta_{s}}{y_{\text{beam}}}\right)\left[\theta(\eta_{0}^{s}-|\eta_{s}|)\right.$
		$\displaystyle\left.+\theta(|\eta_{s}|-\eta_{0}^{s})\exp\left(\frac{\left(|\eta_{s}|-\eta_{0}^{s}\right)^{2}}{2\sigma_{s}^{2}}\right)\right]\,,$		(4)
	$\displaystyle H^{n}_{P/T}(\eta_{s})$	$\displaystyle=\frac{1}{N}\left[\theta(\eta_{s}-|\eta^{n}_{0;P/T}|)\exp\left(-\frac{(\eta_{s}-\eta^{n}_{0;P/T})^{2}}{2\sigma_{n;P/T}^{2}}\right)\right.$
		$\displaystyle+\left.\theta(|\eta^{n}_{0;P/T}|-\eta_{s})\exp\left(-\frac{(\eta_{s}-\eta^{n}_{0;P/T})^{2}}{2\sigma_{n;T/P}^{2}}\right)\right]\,.$		(5)

The parameters in the above envelope functions are fixed by comparing to the final rapidity distributions of charged hadrons. Table 1 shows all the parameters used in the Monte-Carlo Glauber model for initial conditions in Au+Au collisions at different collision energies. Note that in this work, for each collision energy the centrality is classified according to the impact parameter. For each centrality class, we perform 500 event-by-event hydrodynamic simulations. In this work, we set the initial flows in the transverse and space-time rapidity directions to be zeros, i.e., $U_{x}=U_{y}=U_{\eta_{s}}=0$.

To take into account the effect of finite net baryon density at RHIC-BES energies, our integrated event-by-event CLVisc framework includes both energy-momentum and baryon number conservation:

	$\displaystyle\nabla_{\mu}T^{\mu\nu}=0\,,$		(6)
	$\displaystyle\nabla_{\mu}J^{\mu}=0\,,$		(7)

where $\nabla_{\mu}$ represents the covariant derivative operator in the Melin coordinate. The energy-momentum tensor $T^{\mu\nu}$ and net baryon current $J^{\mu}$ take the following form:

	$\displaystyle T^{\mu\nu}=eU^{\mu}U^{\nu}-P\Delta^{\mu\nu}+\pi^{\mu\nu}\,,$		(8)
	$\displaystyle J^{\mu}=nU^{\mu}+V^{\mu}\,,$		(9)

where $e$ is the energy density, $U^{\mu}$ is the flow velocity, $P$ is the pressure, $\Delta^{\mu\nu}=g^{\mu\nu}-U^{\mu}U^{\nu}$, $\pi^{\mu\nu}$ is the shear-stress tensor, $n$ is the net baryon density, and $V^{\mu}$ is baryon diffusion current. In this work, we follow a few recent studies (Shen:2017bsr, ; Akamatsu:2018olk, ; Shen:2020jwv, ; Denicol:2018wdp, ) and neglect the bulk viscosity in the evolution equations for the sake of simplicity and saving the computing time. We would like to point out that the bulk viscosity does have important effects on the medium evolution, especially near the region of the phase transition. It will typically reduce the radial flow of the bulk matter as shown by several previous studies at the LHC energies (Ryu:2015vwa, ; Ryu:2017qzn, ). In order to compensate the effect of bulk viscosity on the radial flow, we have adjusted the initial proper time $\tau_{0}$ to reproduce the final charged particle mean transverse momenta $\langle p_{T}\rangle$. The inclusion of the bulk viscosity in the evolution equations will be left for a future study. The dissipative currents $\pi^{\mu\nu}$ and $V^{\mu}$ are evolved according to the following Israel-Stewart-like equations (Denicol:2018wdp, ):

	$\displaystyle\Delta^{\mu\nu}_{\alpha\beta}D\pi^{\alpha\beta}$	$\displaystyle=-\frac{1}{\tau_{\pi}}\left(\pi^{\mu\nu}-\eta_{v}\sigma^{\mu\nu}\right)$		(10)
		$\displaystyle-\frac{4}{3}\pi^{\mu\nu}\theta-\frac{5}{7}\pi^{\alpha<\mu}\sigma_{\alpha}^{\nu>}+\frac{9}{70}\frac{4}{e+P}\pi^{<\mu}_{\alpha}\pi^{\nu>\alpha}\,,$
	$\displaystyle\Delta^{\mu\nu}DV_{\nu}$	$\displaystyle=-\frac{1}{\tau_{V}}\left(V^{\mu}-\kappa_{B}\bigtriangledown^{\mu}\frac{\mu_{B}}{T}\right)-V^{\mu}\theta-\frac{3}{10}V_{\nu}\sigma^{\mu\nu}\,,$		(11)

where $\theta$ is the expansion rate, $\sigma^{\mu\nu}$ is the symmetric shear tensor, $\eta_{v}$ and $\kappa_{B}$ are the transport coefficients for the evolution of shear tensor and baryon diffusion current. In this work, we consider the specific shear viscosity $C_{\eta_{v}}$ and baryon diffusion coefficient $\kappa_{B}$ as model parameters, which are related to $\eta_{v}$ and $k_{B}$ as follows:

	$\displaystyle C_{\eta_{v}}$	$\displaystyle=\frac{\eta_{v}T}{e+P},$		(12)
	$\displaystyle\kappa_{B}$	$\displaystyle=\frac{C_{B}}{T}n\left(\frac{1}{3}\cot\left(\frac{\mu_{B}}{T}\right)-\frac{nT}{e+P}\right)\,.$		(13)

One can see that the specific shear viscosity $C_{\eta_{v}}$ reduces to $\eta_{v}/s$ for zero net baryon density $n=0$ (used at top RHIC and LHC energies). In this work, we set $C_{\eta_{v}}=0.08$ and $C_{B}=0.4$ as constants for all collision energies. The relaxation times are chosen as $\tau_{\pi}=\frac{5C_{\eta_{v}}}{T}$ and $\tau_{V}=\frac{C_{B}}{T}$.

To solve the hydrodynamic conservation and dissipation equations numerically, the Kurganov-Tadmor (KT) algorithm Pang:2018zzo ; Gale:2013da ; KURGANOV2000241 ; Schenke:2010nt with second-order Runge-Kutta method is utilized to obtain the dynamical evolution of local macroscopic thermodynamic quantities (energy density $e$, net baryon density $n$ and flow velocity $U^{\mu}$). KT algorithm is a finite volume method with the advantages of small numerical viscosity and clear physical interpretation; the change of average conserved quantities in a cell is determined by the flux at the cell interfaces. In order to keep numerical stabilities when solving the dissipative hydrodynamics equations, we compare, at each cell, the maximum components of shear tensor $\pi^{\mu\nu}$ and baryon diffusion current $q^{\mu}$ to the ideal parts of energy-momentum tensor and net baryon current. If ${\rm max}(|\pi^{\mu\nu}|)>T^{\tau\tau}_{\rm ideal}$ or ${\rm max}(|V^{\mu}|)>J^{\tau}_{\rm ideal}$, we locally set $\pi^{\mu\nu}=0$ or $V^{\mu}=0$. Since the cells that need regulations are usually located at the dilute region, such treatment should not affect the major part of medium evolution.

To close the equations of motion for the relativistic hydrodynamics, equation of state of the QCD medium has to be supplied. Several models have implemented the first order phase transition and the critical point in equation of states (Parotto:2018pwx, ). In this work, we use the NEOS-BQS equation of state (Monnai:2019hkn, ; Monnai:2021kgu, ). NEOS-BQS is based on the lattice QCD simulation at high temperature and vanishing net baryon density, and utilizes Taylor expansion method to construct the equation of state at finite net baryon density. At lower energy density, it matches the hadron gas equation of state via a smooth crossover.

When the local energy density drops below the freeze-out energy density (we set $e_{\text{frz}}$= 0.4 GeV/fm³), the Cooper-Frye formula is used to obtain the momentum distribution of particles:

$\displaystyle\frac{dN}{dYp_{T}dp_{T}d\phi}=\frac{g_{i}}{(2\pi)^{3}}\int_{\Sigma}p^{\mu}d\Sigma_{\mu}f_{eq}(1+\delta f_{\pi}+\delta f_{V})\,.$

(14)

In the above equation, $g_{i}$ is the degeneracy for identified hadrons; $d\Sigma_{\mu}$ is the hyper-surface element which is determined from the Cornelius routine Huovinen:2012is ; $f_{\rm eq}$, $\delta f_{\pi}$ and $\delta f_{V}$ are thermal equilibrium distribution and out-of-equilibrium corrections, which take the following forms:

	$\displaystyle f_{\rm eq}$	$\displaystyle=\frac{1}{\exp\left[(p_{\mu}U^{\mu}-B\mu_{B}\right)/T_{f}]\pm 1}\,,$		(15)
	$\displaystyle\delta f_{\pi}(x,p)$	$\displaystyle=(1\pm f^{\text{eq}}(x,p))\frac{p_{\mu}p_{\nu}\pi^{\mu\nu}}{2T^{2}_{f}(e+P)},$		(16)
	$\displaystyle\delta f_{V}(x,p)$	$\displaystyle=(1\pm f^{\text{eq}}(x,p))\left(\frac{n}{e+P}-\frac{B}{U^{\mu}p_{\mu}}\right)\frac{p^{\mu}V_{\mu}}{\kappa_{B}/\tau_{V}},$		(17)

where $T_{f}$ is the chemical freeze-out temperature, $\mu_{B}$ is the net baryon chemical potential, $B$ is the baryon number for the identified baryon, $n$ is the local net baryon density. Note that the above forms of the out-of-equilibrium corrections $\delta f_{\pi}$ and $\delta f_{V}$ are derived from Boltzmann equation via relaxation time approximation (McNelis:2021acu, ). Due to nonzero net baryon density, $T_{f}$ and $\mu_{B}$ are different for each hypersurface cell. This will lead to different particle density when particles are sampled in the comoving frame of the fluid. Note that the thermodynamic variables on the hypersurface are calculated in terms of SMASH hadron resonance gas EOS to be consistent with particle species in SMASH. In this work, we follow a few recent studies (Pang:2018zzo, ; Karpenko:2015xea, ; Shen:2014vra, ; Shen:2014lye, ) and use the Monte-Carlo method to sample the positions and momenta of thermal hadrons according to the Cooper-Frye formula. We use a step function $\theta(f_{\rm eq}+\delta f_{\pi}+\delta f_{V})$ so that the full (equilibrium + non-equilibrium) distribution functions never encounter negative values. Such a treatment should have small effect on the final results since hydrodynamics mainly focuses on the soft physics while large viscous corrections usually appear at large $p_{T}$ region. Recently, the Pratt-Torrieri-Bernhard (PTM) distribution (McNelis:2021acu, ; Pratt:2010jt, ; Du:2019obx, ) and the maximum-entropy distribution (Everett:2021ulz, ) have been proposed to treat the large viscous corrections on the freeze-out hypersurface. In the future, we will implement the PTM distribution or maximum-entropy distribution in our CLVisc framework.

After the particlization of the fluid, the evolution of hadrons is described by the microscopic transport model SMASH (Weil:2016zrk, ; Schafer:2019edr, ; Mohs:2019iee, ; Hammelmann:2019vwd, ; Mohs:2020awg, ) which solves the relativistic Boltzmann equation with elastic collisions, resonance excitations, string excitations and decays for all mesons and baryons up to mass $\sim$ 2 GeV. To improve the statistics, we repeat the sampling and SMASH simulation 2000 times for each single hydrodynamic event, then average over 2000 SMASH events for final analysis.

To reduce the large amount of computing cost for event-by-event hydrodynamic simulations on central processing units (CPUs), CLVisc framework is accelerated on GPUs using OpenCL languages. For example, using the grid size $N_{x}\times N_{y}\times N_{\eta}=201\times 201\times 201$, running a hydrodynamics event with a smooth Glauber initial condition for 0-5% Au+Au collisions at $\sqrt{s_{NN}}=19.6$ GeV on GeForce GTX 3080 GPU typically takes $\sim 1800$ seconds or $\sim 1.6$ seconds per time step. The particle sample module takes $\sim 0.066$ seconds per sample event, and SMASH takes $\sim 2.5$ seconds per afterburner event. Parallelization on GPUs provides a remarkable improvement on computing time for (3+1)-dimensional viscous hydrodynamics simulation, which makes it possible to study the event-by-event properties of QGP at BES energies.

In Figure 1, we show the transverse momentum ($p_{T}$) spectra for identified particles $\pi^{+}$, $K^{+}$, $p$ and $\bar{p}$ in Au+Au collisions at BES energies (7.7-62.4 GeV) with 5 different centrality classes. One can see that our numerical results are in good agreement with STAR data (Adamczyk:2017nof, ) on the $p_{T}$ spectra at mid-rapidity range $|Y|<0.25$ without the contribution from weak decays except at 7.7 GeV. This means that our integrated event-by-event CLVisc hydrodynamics model with finite net baryon density can describe the bulk evolution and radial flow of the QCD matter created in heavy-ion collisions at BES energies.

In order to see the $p_{T}$ spectra and the radial flow effect more clearly, we can look at the mean transverse momenta of identified particles which strongly depend on the slope of the $p_{T}$ spectra. Figure 2 shows the mean $p_{T}$ as a function of the centrality class in Au+Au collisions at BES energies (7.7-62.4 GeV). One can see that our integrated model can describe most of the experimental data from STAR (Adamczyk:2017iwn, ; Abelev:2008ab, ) except the lowest collision energy (7.7 GeV) explored here. For Au+Au collisions at $\sqrt{s_{NN}}=7.7$ GeV, the overlap time is quite large ($\sim 3$ fm) and the pre-equilibrium stage can produce some additional radial flow which is neglected in our model. The use of dynamical initial conditions and pre-equilibrium evolution should improve the model in the future.

In Figure 3, we show the collision energy dependence of the mean $p_{T}$ of identified particles. One can clearly see different radial flow effects for $\pi^{+}$, $K^{+}$ and $p$ (due to different masses): more blue shift effect for more massive particles. As for the collision energy dependence, we find that the mean momenta of $\pi^{+}$, $K^{+}$ and $p$ increase mildly with the collision energy due to larger radial flow.

Anisotropic collective flows, which originate from initial geometric anisotropies and fluctuations, are very important observables for studying the transport properties of QGP produced in relativistic heavy-ion collisions. In this work, we use the $\bm{Q_{n}}$ cumulant method (Bilandzic:2010jr, ) to compute the $n$-th order anisotropic flows. In this method, one first defines the $\bm{Q_{n}}$ vector:

$\displaystyle\bm{Q_{n}}=\sum_{i=1}^{N}e^{in\phi_{i}}\,,$

(18)

where $N$ is the charged multiplicity for each event, and $\phi_{i}$ is the azimuthal angle for the $i$-th particle in the final state. Then the single-event-averaged $2$-particle and $4$-particle correlations can be defined as follows:

	$\displaystyle\left\langle 2\right\rangle$	$\displaystyle=\frac{|Q_{n}|^{2}-N}{N(N-1)}\,,$
	$\displaystyle\left\langle 4\right\rangle$	$\displaystyle=\frac{|Q_{n}|^{4}+|Q_{2n}|^{2}-2\text{Re}[Q_{2n}(Q^{*}_{n})^{2}]}{N(N-1)(N-2)(N-3)}$
		$\displaystyle-2\frac{2(N-2)|Q_{n}|^{2}-N(N-3)}{N(N-1)(N-2)(N-3)}\,.$		(19)

In the above flow correlations $\left\langle 2\right\rangle$ and $\left\langle 4\right\rangle$, the self-correlations and finite particle number effects have been subtracted. Then the 2-particle and 4-particle cumulants can be obtained from $\langle 2\rangle$ and $\langle 4\rangle$ by averaging over many events:

	$\displaystyle C_{n}\{2\}$	$\displaystyle=\left\langle\left\langle 2\right\rangle\right\rangle\,,$		(20)
	$\displaystyle C_{n}\{4\}$	$\displaystyle=\left\langle\left\langle 4\right\rangle\right\rangle-2\left\langle\left\langle 2\right\rangle\right\rangle^{2}\,.$		(21)

Finally, the integrated anisotropic flow coefficients can be calculated as:

	$\displaystyle v_{n}\{2\}$	$\displaystyle=\sqrt{C_{n}\{2\}}\,,$		(22)
	$\displaystyle v_{n}\{4\}$	$\displaystyle=(-C_{n}\{4\})^{\frac{1}{4}}\,.$		(23)

Figure 4 shows the centrality dependence of elliptic flow $v_{2}\{2\}$ and triangle flow $v_{3}\{2\}$ for charged hadrons in Au+Au collisions at 7.7-62.4 GeV using two-particle cumulant method. Here we only consider particles with $p_{T}>0.2$ GeV and $|\eta|<1$ in the flow analysis. One can see that our results are in good agreement with the experimental data from STAR (Adamczyk:2017hdl, ; Adamczyk:2016exq, ). For the elliptic flow, one can clearly see the typical non-monotonic centrality dependence which originates from the combined effect of the elliptic geometry, geometrical fluctuations and the sizes of the collision systems. From central to mid-central to peripheral collisions, the elliptic flow first increases from a finite value in the most central collisions caused by the geometrical fluctuations to the maximum value due to the increased eccentricity, and then decreases due to smaller system sizes. The triangular flow has much weaker dependence on the collision centrality because it mainly originates from the initial state geometrical fluctuations.

In Figure 5, we show elliptic flow $v_{2}\{2\}$ and triangular flow $v_{3}\{2\}$ in most central (0-5%) and mid-central (20-30%) Au+Au collisions as a function of the collision energy (7-62.4 GeV). One can see that for the most central 0-5% collisions, our results on elliptic flow underestimate the experimental data a little bit. One possible reason is the use of constant shear viscosity in our current study. The inclusion of the temperature and baryon chemical potential dependences of the shear viscosity may improve the current result (Shen:2020jwv, ). This will be left for a future study. As for the collision energy dependence, we find that both elliptic and triangular flows increase slightly with the beam energy. This is mainly from the increase of radial flow due to the increase of initial energy density since the eccentricities $\epsilon_{2}$ and $\epsilon_{3}$ have very weak dependence on collision energy as shown in the figure.

The fluctuations of anisotropic flows are also valuable tools to probe the initial state fluctuations and the transport properties of the QGP. In Figure 6, we show the multi-particle cumulant ratio $v_{2}\{4\}/v_{2}\{2\}$ as a function of centrality for different collision energies (7-62.4 GeV). The deviation of this cumulant ratio from unity quantifies the relative flow fluctuations. One can clearly see the non-monotonic dependence of the cumulant ratio $v_{2}\{4\}/v_{2}\{2\}$ on the collision centrality. As one goes from central to mid-central to peripheral collisions, the multi-particle cumulant ratio $v_{2}\{4\}/v_{2}\{2\}$ first increases, and then decreases. This means that the relative fluctuations are smallest in mid-central collisions, and are larger in central and peripheral collisions. This is because the initial collision geometry dominates elliptic flow in mid-central collisions, while the fluctuations dominate elliptic flow in central and peripheral collisions. Another interesting observation is that the multi-particle cumulant ratio $v_{2}\{4\}/v_{2}\{2\}$ has weak collision energy dependence, consistent with the preliminary STAR data (Magdy:2018itt, ; Niseemtalk, ).