Effects of bulk viscosity and hadronic rescattering in heavy ion collisions at RHIC and LHC

TL;DR

This paper uses a hybrid model to study how bulk viscosity and hadronic rescattering affect observable outcomes in ultra-relativistic heavy ion collisions at RHIC and LHC, highlighting the importance of bulk viscosity.

Contribution

It introduces a hybrid modeling approach combining IP-Glasma, viscous hydrodynamics, and microscopic transport to analyze the effects of bulk viscosity and rescattering.

Findings

Bulk viscosity significantly influences radial flow and azimuthal anisotropy.

The model reproduces experimental data on flow observables.

Dependence of results on the switching temperature to microscopic transport.

Abstract

We describe ultra-relativistic heavy ion collisions at RHIC and the LHC with a hybrid model using the IP-Glasma model for the earliest stage and viscous hydrodynamics and microscopic transport for the later stages of the collision. We demonstrate that within this framework the bulk viscosity of the plasma plays an important role in describing the experimentally observed radial flow and azimuthal anisotropy simultaneously. We further investigate the dependence of observables on the temperature below which we employ the microscopic transport description.

Equations

$$⟨{\rho }_i^a({\text{x}}_T^{\prime }){\rho }_j^b({\text{x}}_T^{\prime \prime })⟩=g^2{\mu }_A^2({\text{x}}_T^{\prime }){\delta }^{ab}\frac{{\delta }^{ij}}{N_L}{\delta }^{(2)}({\text{x}}_T^{\prime }-{\text{x}}_T^{\prime \prime }),$$

$$Q_{s,A}^2({\text{x}}_T)$$

$${\mu }^2(r^2)$$

$$x{f}_g(x,{\mu }^2={\mu }_0^2)=A_g{x}^{-{\lambda }_g}(1-x{)}^{5.6}.$$

$$A_{(1,2)}^i=-\frac{i}{g}{U}_{(1,2)}{\partial }_i{U}_{(1,2)}^{†}({\text{x}}_T)$$

$$U_{(1,2)}({\text{x}}_T)=\mathcal{P}\exp \left[-ig\int d{x}^{(+,-)}\frac{{\rho }_{(1,2)}({\text{x}}_T,x^{(+,-)})}{{\nabla }_T^2-m^2}\right].$$

$$A^i(\tau \to 0^{+})$$

$$A^{\eta }(\tau \to 0^{+})$$

$${\partial }_{\mu }{F}^{\mu \nu }-ig[A_{\mu },F^{\mu \nu }]=0$$

$$F^{\mu \nu }={\partial }^{\mu }{A}^{\nu }-{\partial }^{\nu }{A}^{\mu }-ig[A^{\mu },A^{\nu }].$$

$$T^{\mu \nu }=-2\text{Tr}\left(F_{\alpha }^{\mu }{F}^{\nu \alpha }\right)+\frac{1}{2}{g}^{\mu \nu }\text{Tr}\left(F^{\alpha \beta }{F}_{\alpha \beta }\right)$$

$$T_{\nu }^{\mu }{u}^{\nu }=ϵu^{\mu }$$

$${\partial }_{\mu }{T}^{\mu \nu }=0.$$

$$T^{\mu \nu }=ϵu^{\mu }{u}^{\nu }-\left(P+\Pi \right){\Delta }^{\mu \nu }+{\pi }^{\mu \nu }$$

$${\tau }_{\Pi }\dot{\Pi }+\Pi$$

$${\tau }_{\pi }{\dot{\pi }}^{⟨\mu \nu ⟩}+{\pi }^{\mu \nu }$$

$${\tau }_{\pi }$$

$${\tau }_{\Pi }$$

$$\frac{{\delta }_{\pi \pi }}{{\tau }_{\pi }}$$

$$\frac{{\tau }_{\pi \pi }}{{\tau }_{\pi }}$$

$$\frac{{\lambda }_{\pi \Pi }}{{\tau }_{\pi }}$$

$$\frac{{\delta }_{\Pi \Pi }}{{\tau }_{\Pi }}$$

$$\frac{{\lambda }_{\Pi \pi }}{{\tau }_{\Pi }}$$

$${\phi }_7$$

$$\frac{dN}{d^3\text{p}}$$

$$\begin{array}[]{l}\displaystyle\bar{N}|_{\scriptsize\textrm{1-cell}}\vspace{2pt}\\ \displaystyle=\left\{\begin{array}[]{ll}\displaystyle\left[n_{0}(x)+\delta n_{\scriptsize\textrm{bulk}}(x)\right]u^{\mu}\Delta\Sigma_{\mu}&\quad\textrm{if}\quad u^{\mu}\Delta\Sigma_{\mu}\geq 0\vspace{2pt}\\ 0&\quad\textrm{otherwise}\end{array}\right.\end{array}$$

$$n_0(x)$$

$$\delta n_{\text{bulk}}(x)$$

$$⟨N^2⟩-⟨N{⟩}^2$$

$$\ne$$

$${\frac{dN}{d^3\text{p}}|}_{\text{1-cell}}=\frac{d}{(2\pi {)}^3}\left[f_0+\delta f_{\text{shear}}+\delta f_{\text{bulk}}\right]\frac{p^{\mu }\Delta {\Sigma }_{\mu }}{E_{\text{p}}}$$

$$f_0$$

$$\delta f_{\text{shear}}$$

$$\delta f_{\text{bulk}}$$

$$\frac{1}{C_{\text{bulk}}}$$

$$f(P,{\pi }^{\mu \nu },\Pi )$$

Full text

I Introduction

Ultra-relativistic heavy ion collisions carried out at the Relativistic Heavy-ion Collider (RHIC) and the large Hadron Collider (LHC) are unequaled tools to study the many-body properties of quantum chromodynamics (QCD), in particular its high-temperature deconfined phase known as the quark-gluon plasma (QGP) gale:2013 ; deSouza:2015ena . Since the QGP is only produced for a very short time and cannot be observed directly, extracting its properties from heavy ion measurements is a major challenge that requires modelling the many stages of the collision: the pre-equilibrium dynamics of the system, the rapid expansion and cooling of the QGP and the dynamics of the dilute hadronic matter that is eventually measured by the experiments.

In the past decade, hydrodynamic models have been applied with great success to describe the distribution of soft hadrons produced in heavy ion collisions at RHIC and the LHC gale:2013 ; deSouza:2015ena . The foremost experimental discovery made using these models was that the QGP displays remarkable transport properties, with one of the smallest shear viscosity to entropy density ratio ever observed Gyulassy:2004zy ; whitepaper1 ; whitepaper2 ; whitepaper3 ; whitepaper4 . A more precise determination of the transport properties of QCD matter, including their non-trivial temperature dependence, is one of the primary goals of the heavy ion research program.

For a long time, shear viscosity was considered to be the dominant source of dissipation for the QGP produced in heavy ion collisions Romatschke:2007mq ; MUSIC ; Song:2010aq ; Song:2011qa ; Arnold:2006bv . Nevertheless, there are theoretical indications that bulk viscosity can become large around the QCD crossover region Meyer:2007dy ; Karsch:2007jc ; Buchel:2007mf ; Kharzeev:2008tra ; NoronhaHostler:2008ju and can significantly affect the evolution of the QGP Torrieri:2008ip ; Rajagopal:2009yw ; Denicol:2009am ; Habich:2014tpa ; Denicol:2015bpa . Early investigations on the effect of bulk viscosity using realistic hydrodynamic simulations often assumed small values for this transport coefficient Bozek:2009dw ; Bozek:2012qs ; Song:2009rh and found modest effects. Other studies focused on the effects of dissipative corrections due to bulk viscosity in the particlization Huovinen:2012is of the hadron resonance gas Monnai:2009ad ; Dusling:2011fd ; Noronha-Hostler:2013gga ; Noronha-Hostler:2014dqa . Whether a large bulk viscosity can be reconciled with the current theoretical description of heavy ion collisions is a topic of great interest to the field.

Recent calculations done in Ref. Ryu:2015vwa have addressed this issue within a modern hydrodynamic description, finding that a large bulk viscosity around the phase transition region is essential to describe simultaneously the multiplicity and average transverse momentum of charged hadrons. This finding was made using IP-Glasma initial conditions Schenke:2012wb , second order hydrodynamic equations MUSIC ; Denicol:2014vaa , and a transport description of the late stages of the collision UrQMD1 ; Bleicher:1999xi . Similar conclusions about the importance of bulk viscosity have since been reached by calculations employing different initial state models that, similarly to IP-Glasma, also exhibit large sub-nucleonic energy density fluctuations Denicol:2015nhu ; Bernhard:2016tnd .

The goal of this paper is to expand on the results presented in Ref. Ryu:2015vwa , offering a more detailed overview on the effect of bulk viscosity on other heavy ion observables and different collision energies. In particular, the effects of late stage hadronic rescattering Bass:2000hy ; Teaney:2001fl will be discussed in greater details. These have been investigated in a number of previous publications Song:2010aq ; Song:2011qa ; Song:2013qma ; Zhu:2015dfa ; Hirano:2005xf ; Hirano:2007ei ; Hirano:2010jg ; Takeuchi:2015ana ; Nonaka:2007hy ; Petersen:2008in ; Nonaka:2010hc ; Petersen:2010tr , where they were found to be especially important to provide a reasonable description of the hadronic chemistry of heavy ion collisions, specially for heavier baryons. Here, we shall perform a systematic study of the effects of the switching temperature between the hydrodynamic simulation and the transport model, showing that this parameter has a significant effect on the momentum distribution of protons and multi-strange hadrons.

The rest of this paper is organized as follows. In Section II, we explain each of the components of our model and show how they are combined to give an integrated description. In Section III, we compare the results of our calculations with experimental data, focusing on the effect of bulk viscosity and hadronic rescattering. We show that a finite bulk viscosity resolves the tension between the observed multiplicity and the average transverse momentum. We also show that the hadronic cascade is an important ingredient for the description of hadronic chemistry. We summarize our results and discuss implications in Section IV.

II Model

The theoretical framework used in the present paper can be divided in four parts: the pre-equilibrium dynamics described with the IP-Glasma model, the hydrodynamical evolution, the transition from fluid to particles and the final hadronic transport.

II.1 Pre-equilibrium with IP-Glasma

The IP-Glasma model Schenke:2012wb describes the pre-equilibrium dynamics of a large number of low-$x$ gluons by the classical Yang-Mills equation. The high-$x$ partons serve as the color sources for the initial gluon fields before the collision. For each nucleus, color charges $\rho_{i}^{a}(\textbf{x}^{\prime}_{T})$ with color index $a$ at a lattice site $i$ in the $x^{\pm}$-direction are sampled according to a Gaussian distribution, satisfying

[TABLE]

where $N_{L}$ is the number of lattice sites in the $x^{\pm}$-direction. The value $N_{L}=100$ is used in this work. The average color charge density per unit transverse area $g^{2}\mu_{A}^{2}(\textbf{x}_{T})$ is proportional to the saturation scale $Q_{s,A}^{2}(\textbf{x}_{T})$ determined in the impact parameter dependent saturation model (IP-Sat) Bartels:2002msat ; Kowalski:2003ipsat

[TABLE]

where $f_{g}(x,\mu^{2})$ is the gluon distribution function in a nucleon and $\mu_{0}$ is a momentum scale at which the gluon distribution has a form of Bartels:2002msat

[TABLE]

The gluon distribution at an arbitrary momentum scale $\mu>\mu_{0}$ is obtained from DGLAP evolution. Each nucleon is assumed to have a Gaussian shape with width $\sigma_{0}$ in the transverse plane. The parameters are determined to fit HERA deep inelastic scattering data Kowalski:2006dip . In a $q\bar{q}+p$ scattering, $r$ is the size of the $q\bar{q}$ dipole and corresponds to the spatial scale of the probe.

The positions of nucleons inside a nucleus are sampled according to the Wood-Saxon distribution and $\textbf{x}_{T,i}$ is the position of the $i$-th nucleon in the transverse plane. Once the color charge distribution is determined, we solve the classical Yang-Mills equation to obtain the gluon field in each nucleus

[TABLE]

where the Wilson line $U$ is

[TABLE]

The subscript $(1)$ and $(2)$ indicates projectile and target quantities, respectively. Since the fluctuation scale of $\rho_{(1,2)}$ is $\sim Q_{s}$, so is the fluctuation scale in $A^{i}_{(1,2)}$.

The gluon field right after the collision ($\tau\to 0^{+}$) is given by Krasnitz:1999wc ; Krasnitz:2000gz ; Krasnitz:2002mn

[TABLE]

For $\tau>0$, we evolve the gluon field according to the Yang-Mills equation

[TABLE]

where the field strength tensor is given as usual by

[TABLE]

After evolving the gluon field up to $\tau_{0}=0.4\,\textrm{fm}$, the energy-momentum tensor is formed out of the field strength tensor

[TABLE]

where the trace is over color in the fundamental representation. The time-like eigenvalue of $T^{\mu}_{\ \nu}=T^{\mu\lambda}g_{\lambda\nu}$ provides the local energy density and the flow velocity

[TABLE]

of IP-Glasma at $\tau_{0}$.

The normalization of the energy-momentum tensor in IP-Glasma is not fully constrained, owing to freedom in the choice of $\alpha_{s}$, as noted in Ref. Schenke:2012fw . This normalization can be fixed by comparing the results of the hydrodynamical simulation with charged hadron multiplicity measurements. This is the procedure adopted in this work. Effectively, this translates into a normalization of the energy density $\epsilon$ of IP-Glasma. The initial flow $u^{\mu}$ is unaffected by this normalization of $T^{\mu\nu}$. Note that the shear stress tensor of IP-Glasma is not currently used to initialize the hydrodynamic simulation, where $\pi^{\mu\nu}$ is initialized to zero.

II.2 Second-order relativistic hydrodynamics

with shear and bulk viscosity

The main hydrodynamic equations are the conservation laws of net-charge, energy, and momentum. Since we only aim to describe the matter produced in the mid-rapidity region at high collision energies, the net-charge can be approximated to be zero and we are only required to solve the continuity equation for $T^{\mu\nu}$,

[TABLE]

In a viscous fluid, the energy-momentum tensor $T^{\mu\nu}$ is decomposed in terms of the velocity field as

[TABLE]

where $\Delta^{\mu\nu}\equiv g^{\mu\nu}-u^{\mu}u^{\nu}$, $P$ is the thermodynamic pressure, $\Pi$ is the bulk viscous pressure, and $\pi^{\mu\nu}$ is the shear stress tensor. The relation between $\epsilon$ and $P$ is given by an equation of state, $P(\epsilon)$. In this work we use the equation of state constructed from a hadronic resonance gas and lattice calculation Huovinen:2009yb .

The time evolution of the bulk and shear viscous corrections, driven by the expansion rate $\theta=\nabla_{\mu}u^{\mu}$ and the shear tensor $\sigma^{\mu\nu}=\frac{1}{2}\left[\nabla^{\mu}u^{\nu}+\nabla^{\nu}u^{\mu}-\frac{2}{3}\Delta^{\mu\nu}(\nabla_{\alpha}u^{\alpha})\right]$, in which $\nabla_{\mu}=(g_{\mu\nu}-u_{\mu}u_{\nu})\partial^{\nu}$, are given by the equations

[TABLE]

which follows from the 14-moment approximation of the Boltzmann equation  Denicol:2012cn ; Denicol:2014vaa . The first-order transport coefficients $\eta$ and $\zeta$ are the shear and bulk viscosities, respectively. The shear and bulk relaxation time, $\tau_{\pi}$ and $\tau_{\Pi}$ are found to be

[TABLE]

where $c_{s}$ is the speed of sound. The second-order transport coefficients are related to the relaxation time through the relations

[TABLE]

These values were first obtained in Denicol:2014vaa . Once we have $\eta$ and $\zeta$ as functions of temperature, it is possible to find the temperature dependence of the relaxation times and the remaining second-order transport coefficients. The shear viscosity over the entropy density ratio $\eta/s$ is set to be constant in this work. Data-driven determination of the temperature-dependent $\eta/s$ via Bayesian analysis is performed in Bernhard:2016tnd . The temperature dependence of the bulk viscosity to entropy density ratio $\zeta/s$ is fixed as shown in Fig. 1 and the temperature where bulk viscosity peaks is set to be $T_{\scriptsize\textrm{peak}}=180\,\textrm{MeV}$ based on the transition temperature in the equation of state.

II.3 Transition from hydrodynamics to transport theory

During the hydrodynamic evolution, the system becomes gradually more dilute and, at some point, a hydrodynamic description will break down. Nevertheless the system is still interacting and the subsequent dynamics must be described in another framework, transport theory for example. Using transport generally means that fluid elements must be converted into hadronic degrees of freedom, which will then be described using a hadronic kinetic theory simulation. In principle, this matching between degrees of freedom must be performed in a space-time region in which both hydrodynamics and kinetic theory are within their domain of applicability.

One way to transition from hydrodynamics to transport theory is to do so when the expansion rate of the fluid is moderately (but not significantly) smaller than the mean-free path of the hadrons composing the fluid. In this way, the system is interacting enough for hydrodynamics to apply and dilute enough for the Boltzmann equation to be applicable. In practice, implementing such procedure can be rather complicated, since it requires extensive knowledge about the interactions among the hadrons. In this work, we adopt the common simplification of approximating the switching hypersurface as a constant temperature hypersurface, with the switching temperature $T_{\rm sw}$ becoming one of the free parameters of our model. As we shall discuss in the following sections, this parameter will be determined by optimizing the fit of identified charged hadron multiplicity, especially protons.

It is important to emphasize that $T_{\rm sw}$ is an effective parameter of the model, that is supposed to describe the more complicated physics of the transition from the description of dense to dilute systems. In this sense, there is no reason to expect such temperature to remain the same as one changes the collision energy (from RHIC to LHC) or even centrality class. One of our conclusions will be that the switching temperature at RHIC (=165 MeV) is larger than the switching temperature at the LHC (=145 MeV). This difference may be due to the fact that systems produced at LHC energies have more entropy and, consequently, are more long-lived than the ones produced at RHIC. Naturally, a more precise explanation for the switching parameters we extract can only be obtained by improving the model, taking into account a more realistic transitioning to the transport phase with improved viscous correction to the distribution function.

We now present the details of how we switch from hydrodynamics to hadronic transport. On isothermal hypersurfaces with constant switching temperature $T_{\rm sw}$, we sample particles with degeneracy $d$ and mass $m$ according to the Cooper-Frye formula Cooper:1974mv

[TABLE]

where $E_{\scriptsize\textbf{p}}$ satisfies $E_{\scriptsize\textbf{p}}^{2}=\textbf{p}^{2}+m^{2}$. The normal vector $d^{3}\Sigma_{\mu}$ is an exterior product of three displacement vectors tangential to the hypersurface. In our simulation, we construct the hypersurface from tetrahedra MUSIC and sample hadrons at each of the grid locations $x$. We approximate the probability distribution of the number of particles to be a Poisson distribution whose average value is given by

[TABLE]

where the number density at thermal equilibrium $n_{0}$ and the bulk viscous correction $\delta n_{\scriptsize\textrm{bulk}}$ are given by

[TABLE]

The shear tensor correction does not induce a change in the number density because of its spin-2 structure (see Eq.(32) below), which is orthogonal to any scalar.

It is understood that quantum thermal distributions are not Poissonian, since

[TABLE]

Nevertheless, within the range of switching temperature considered in this work ($135\,\hbox{MeV}\leq T_{\rm sw}\leq 165\,\hbox{MeV}$), we verified that these quantum effects are less than $10\,\%$ for pions and less than $1\,\%$ for heavier hadrons. Using the Poisson distribution is therefore a reasonable approximation.

After we determine the number of particles in each cell, we sample the momentum of each particle according to the following prescription Huovinen:2012is

[TABLE]

if $(f_{0}+\delta f_{\scriptsize\textrm{shear}}+\delta f_{\scriptsize\textrm{bulk}})>0$ and $p^{\mu}\Delta\Sigma_{\mu}>0$. Otherwise, $\left.dN/d^{3}\textbf{p}\right|_{\scriptsize\textrm{1-cell}}=0$. In our simulation, the deviation from particle spectra given by Riemann integration of Eq.(24) is less than 5% for $p_{T}<2\,\textrm{GeV}$. Therefore, as long as the soft physics is concerned, Eq.(30) is an adequate implementation of the Cooper-Frye transition.

The explicit expressions for the equilibrium distribution functions and the shear Dusling:2010rel and bulk Bozek:2009dw ; Paquet:2015lta viscous corrections are

[TABLE]

where

[TABLE]

and the flow velocity $u^{\mu}$ and temperature $T$ on the hypersurface are determined from the hydrodynamic evolution. The summation is over hadronic species.

II.4 Microscopic transport UrQMD as afterburner

The sampled particles are propagated in UrQMD (version 3.4) UrQMD1 ; Bleicher:1999xi , which simulates interactions of hadrons and resonances with masses up to $2.25\,\textrm{GeV}$. These interactions include inelastic processes through resonance scattering, $B\bar{B}$ annihilation and string excitation, as well as elastic scatterings. Whenever experimental data are available, the hadronic cross sections in UrQMD are based on the data. When measurements are not available, cross sections are extrapolated from other processes based on detailed balance and the additive quark model. Using UrQMD as afterburner allows for a more realistic description of the late stage of the collision, where the mean free path is not short compared to the macroscopic scale given by system size or expansion rate.

Note that while baryon-antibaryon annihilation is included in UrQMD, pair creation is not. This is because $B\bar{B}$ predominantly annihilate into multiple pions, but the opposite channel, which would involve the simultaneous interaction of multiple hadrons, is not currently supported by UrQMD. Consequently, all $B\bar{B}$ pairs in the system originate from the Cooper-Frye procedure and, strictly speaking, detailed balance is not obeyed. This violation of detailed balance is not expected to be a major issue: previous works such as Ref. Pan:2014caa have shown that the contribution of baryon-antibaryon creation is considerably smaller than that of $B\bar{B}$ annihilation. It is also possible to make the argument that since the system is expanding and the mean free path is comparable to the macroscopic scale, there are more baryons and anti-baryons than there would be in local thermal equilibrium. This is because the system does not have enough interactions to reach equilibrium. The excess of mesons over thermal equilibrium is less significant owing to the lower masses. Therefore, one can expect that $B\overline{B}$ annihilation will be more frequent than the inverse process in the evolution toward equilibrium. In this sense, switching from hydrodynamics to transport coincides with the point where the $B\bar{B}$ annihilation becomes dominant.

III Results and discussion

In this section, we discuss the results of our simulations for Au-Au collisions at RHIC ($\sqrt{s_{NN}}=200$ GeV) and Pb-Pb collisions at the LHC ($\sqrt{s_{NN}}=2.76$ TeV) where a wide set of measurements are available. The centrality classes $0-5\%$, $10-20\%$, $20-30\%$ and $30-40\%$ are considered. We highlight the effect of bulk viscosity and the importance of the hadronic rescattering stage. We show that our approach is capable of describing a large number of hadronic observables consistently with a fixed set of parameters. The main parameters in this work are the switching temperature between the hydrodynamic expansion and the afterburner UrQMD, $T_{\rm sw}$, and the value of the effective shear viscosity over the entropy density ratio $\eta/s$.

III.1 Integrated observables

Observables integrated over the transverse momentum $p_{T}$ generally have a reduced sensitivity to out-of-equilibrium corrections of the hadronic momentum distribution ($\delta f_{\scriptsize\textrm{shear}}$) because of shear viscosity, compared to $p_{T}$-differential observables. Therefore, the multiplicity $\left.dN/dy\right|_{|y|<0.5}$, the mean transverse momentum $\langle p_{T}\rangle$, and the $p_{T}$-integrated anisotropic flow coefficients $v_{n}$ are investigated first. The anisotropic flow $v_{n}$ is computed using the multi-particle cumulant method based on the flow correlations among particles as in Ref. Bilandzic:2010jr .

Figure 2 shows the multiplicity and average $p_{T}$ for pions, kaons and protons, as well as the charged hadron anisotropic flow coefficients $v_{2,3,4}$, from central to semi-peripheral centrality bins. The charged hadron $v_{2}$ is shown for the two- and four- particle cumulants. Calculations that include both shear and bulk viscosities are in the left panels, while the calculations with only shear viscosity are presented in the right panels. In both cases, the value of $\eta/s$ was adjusted such that the measured charged hadron $v_{n}$ is reproduced. A value of $\eta/s=0.095$ is used when both bulk and shear viscosity are present, while a larger value of $\eta/s=0.16$ is necessary in absence of bulk viscosity. This important effect of bulk viscosity on phenomenological extractions of $\eta/s$ had been quantified previously in Ref. Ryu:2015vwa .

Another significant effect of bulk viscosity is a considerable suppression of the average transverse momentum of hadrons, which can be seen by comparing the left and right hand sides of Fig. 2. The change in average $p_{T}$ is significant for all hadron species (pions, kaons, protons). Bulk viscosity is essential for a simultaneous description of the multiplicity and $\langle p_{T}\rangle$ of hadrons when IP-Glasma initial conditions are used: without bulk viscosity, the system expands too rapidly, leading to a larger hydrodynamic transverse flow than suggested by average $p_{T}$ measurements. Bulk viscosity improves the agreement with data by acting as a resistance to expansion, reducing the transverse flow of the system. Besides this change in the plasma expansion, part of the modification of $\langle p_{T}\rangle$ is from the effect of bulk viscosity on the hadronic momentum distribution — $\delta f_{\scriptsize\textrm{bulk}}$ given by Eq.(33). The average $p_{T}$ is actually decreased by $\delta f_{\scriptsize\textrm{bulk}}$. If $\delta f_{\scriptsize\textrm{bulk}}$ were smaller, a similar suppression in $\langle p_{T}\rangle$ could be achieved with a larger bulk viscosity. More details about the effect of $\delta f_{\scriptsize\textrm{bulk}}$ on integrated hadronic observables are presented in Appendix A. We highlight here that $\delta f_{\scriptsize\textrm{bulk}}$ has a small effect on the $v_{n}$ of charged hadrons and the multiplicity of pions and kaons.

The switching temperature $T_{\rm sw}$ between the hydrodynamic simulation and the hadronic afterburner is $145$ MeV for the calculations presented in Fig. 2. This choice can be understood from Fig. 3 (solid lines), which shows the dependence on $T_{\rm sw}$ of the multiplicity and average transverse momentum of identified hadrons (the charged hadron momentum anisotropies have a small dependence on $T_{\rm sw}$ — this is discussed in more details below). The multiplicities of pions and kaons are shown to have a weak dependence on $T_{\rm sw}$, while protons are much more sensitive to this parameter. The value of $T_{\rm sw}$ around $145$ MeV leads to the best agreement with ALICE measurements for the proton multiplicity.

Figure 3 (solid lines) also shows that dependence of the average transverse momentum on $T_{\rm sw}$ is mild for all three hadrons species. We verified that a similar $T_{\scriptsize\textrm{sw}}$ dependence was found for more peripheral collisions, up to the $30-40\%$ centrality class, for both the multiplicity and $\langle p_{T}\rangle$. We further verified that the dependence on the switching temperature $T_{\scriptsize\textrm{sw}}$ of these same observables is very similar with and without bulk viscosity. The effect of $T_{\scriptsize\textrm{sw}}$ on the identified hadron $\langle p_{T}\rangle$ is thus small compared to the effect of bulk viscosity. In consequence, we emphasize that in absence of bulk viscosity, it would not be possible to obtain a good agreement with identified hadron $\langle p_{T}\rangle$ by changing the value of $T_{\scriptsize\textrm{sw}}$.

Also shown in Fig. 3 is how hadronic rescattering affects the $T_{\scriptsize\textrm{sw}}$ dependence of the multiplicity and average transverse momentum. It is found that calculations that include hadronic decays but not hadronic interactions (dashed lines) have a larger slope in $T_{\scriptsize\textrm{sw}}$ than those that include both hadronic decays and rescattering (solid lines). This means that hadronic rescattering reduces the dependence on $T_{\scriptsize\textrm{sw}}$. The effect is fairly small for pions and kaons, but significant for protons. Heavier baryons, shown in Fig. 4, have a similar dependence on $T_{\scriptsize\textrm{sw}}$ as protons. Since hydrodynamics describes an interacting medium, it is indeed expected that the transition between hydrodynamics and UrQMD will be smoother — if not necessarily smooth — when hadronic rescattering is included. The larger dependence of protons and heavier hadrons on $T_{\scriptsize\textrm{sw}}$ can be seen as a systematic uncertainty of our model for observables involving these hadrons.

The effect of hadronic rescattering on our multiplicity and mean $p_{T}$ calculations is shown again in Fig. 5, this time as a function of the centrality class, and compared against ALICE data. The switching temperature is fixed $T_{\scriptsize\textrm{sw}}=145$ MeV, which as explained above provides a good description of the proton multiplicity. It can be seen in Fig. 5 that the effect of hadronic rescattering is very similar across centralities. We also show in Fig. 5 the explicit effect of a subset of hadronic rescatterings, namely baryon-antibaryon annihilation. What is interesting about $B\bar{B}$ annihilations is that they represent most of the change in the multiplicity of protons due to hadronic rescattering, although they have a very small effect on the average transverse momentum of protons. This result is in general agreement with the observations made in Song:2013qma .

In Fig. 6, we show the results of our calculations at RHIC by comparing with STAR measurements Adams:2004bi ; Abelev:2008ab ; Adamczyk:2013waa . Like it was done with LHC calculations, the value of $T_{\rm sw}$ was adjusted so as to provide a good description of the proton multiplicity, while the value of $\eta/s$ was fixed using the charged hadron $v_{2}$. The values of $T_{\rm sw}=165$ MeV and $\eta/s=0.06$ were found to provide good agreement with the respective measurements. The quality of agreement with measurements can be seen to be similar to that found at the LHC (c.f. Fig. 2). We verified that the dependence on $T_{\scriptsize\textrm{sw}}$ of the multiplicity and average transverse momentum of identified hadrons is very similar at RHIC as at the LHC.

To conclude this section on integrated observables, in Fig. 7 we investigate the effect of the switching temperature $T_{\scriptsize\textrm{sw}}$ between hydrodynamics and UrQMD on the momentum anisotropy $v_{2}\{2\}$ of integrated charged hadrons, at RHIC as well as at the LHC. The upper curves correspond to peripheral 30-40% collisions, and the lower curves to central 0-5% collisions. The solid line corresponding to the calculation with hadronic decays and rescattering, and the dashed line including only hadronic decays but not rescattering. As observed previously for the multiplicity and average transverse momentum, the inclusion of hadronic rescattering reduces significantly the observable’s dependence on $T_{\scriptsize\textrm{sw}}$ (i.e. solid lines flatter than dashed ones). At both RHIC and LHC, hadronic rescattering increases $v_{n}$, which is consistent with the effect of rescattering observed on pions in Ref. Hirano:2007ei . This increase is larger at RHIC than at the LHC, and is also larger for peripheral events than for central ones. Our understanding is that this is a consequence of the different lifetime of the hadronic transport phase compared to the hydrodynamic expansion for the different centralities and collision energies, as well as a consequence of how isotropic each system is at the transition between hydrodynamics and transport.

III.2 Differential observables

In this section, we examine $p_{T}$ differential observables. At this point, all model parameters have already been fixed with integrated observables both at LHC and RHIC energies.

The $p_{T}$ differential $v_{n}\{2\}$ of $n=2,3,4$ of charged hadrons are compared with the ALICE ALICE:2011ab and the CMS CMS:2013v2 ; CMS:2014vn data in Fig. 8. Note that the $p_{T}$-differential $v_{n}$ is evaluated from the azimuthal correlation between particles of interest and reference flow particles, given that the particles of interest are those in specific $p_{T}$ bins Bilandzic:2010jr . Although $v_{2}\{2\}$ deviates from the data at high $p_{T}$, especially when compared with the ALICE measurements, our calculation shows a reasonable agreement with data for $p_{T}\lesssim 1$ GeV, where we have the most particles.

We next turn to identified hadron observables at LHC energies. The $p_{T}$-differential spectra of pions, kaons and protons are shown in Fig. 9, with (solid line) and without (dashed line) the effect of hadronic rescattering, for four different centralities. Calculations that include hadronic rescattering agree very well with measurements for the most central collisions ($0-5\%$), for all three hadron species. Tension with data appears and increases in more peripheral centralities, especially in kaons and protons, but also in pions at $p_{T}$ above $1.5-2$ GeV. As expected from the discussion of integrated observables, the hadronic transport phase has a minor effect on the pion spectra, which is slightly hardened at $p_{T}>2$ GeV. The kaon spectra get flatter resulting in a better agreement with the experimental measurement. A more significant effect of rescattering is seen in the proton spectra: the low $p_{T}$ parts of the spectra is reduced in the transport phase owing to $B\bar{B}$ annihilations while hadronic rescattering shifts more protons to higher $p_{T}$. This shows once again that the inclusion of the hadronic transport phase is important to describe the measured proton spectra at the LHC.

Figure 10 shows identified particle elliptic flow coefficients at the LHC, with measurements from the ALICE collaboration ALICE:2014v2id . Comparing the simulation results with (solid lines) and without (dashed lines) hadronic rescattering, we find once more that pions and kaons $v_{2}(p_{T})$ are largely insensitive to rescattering. On the other hand, hadronic rescattering has a large effect on the proton $v_{2}(p_{T})$, which is considerably decreased by hadronic interactions. Even though $v_{2}$ around the mean $p_{T}$ is well reproduced, our calculations overestimate the $v_{2}$ of pions and kaons at higher $p_{T}$. We highlight that tension with ALICE measurements was also observed at high $p_{T}$ for the $v_{2}$ of charged hadrons shown in Fig. 8. We note that tension with measurements at high $p_{T}$ is less worrying than in lower regions of transverse momenta, since this region of $p_{T}$ is more sensitive to uncertainties in the viscous corrections to the hadron distribution function ($\delta f$), as well as potential contribution from recombination with (mini-)jet shower partons. Nevertheless, there still seems to be room for improvement at lower $p_{T}$ in our description of identified hadron $v_{n}$.

The $p_{T}$ spectra and $v_{2}$ of strange baryons are shown in Fig.  11 and compared with the ALICE data ALICE:2014mstr ; Abelev:2013Ks0L ; ALICE:2014v2id . The $p_{T}$ dependence of the spectra of $\Lambda$, $\Xi$ and $\Omega$ is described well, although deviations of up to $20\%$ are observed in the normalization. The effect of hadronic rescattering, which suppresses the $p_{T}$ spectra more at low $p_{T}$, is consistent with the decrease in multiplicity and the increase in average transverse momentum seen in the previous section. We consider the level of agreement with experimental data to be acceptable considering the non-negligible dependence of heavier hadrons on the switching temperature between hydrodynamics and UrQMD shown previously in Fig. 3.

For all three heavy strange baryons, our calculation overestimates the $v_{2}(p_{T})$. Previous studies Zhu:2015dfa ; Takeuchi:2015ana , also based on a hybrid approach with isothermal particlization, found some tension with hyperons as well, although we highlight that comparisons with these previous models is not straightforward because of differences in the hydrodynamic modeling (e.g. initial conditions). Once again, since heavy hadrons have been shown in this work to be especially sensitive to the transition between hydrodynamics and the afterburner, this tension is not unexpected. There have been proposals in the literature that strange hadrons may chemically freeze out earlier than non-strange particles vanHecke:1999jh ; Arbex:2001vx ; Chatterjee:2013yga . We cannot necessarily conclude this from our investigations, but we can say that improvements in the transition between hydrodynamics and hadronic transport are important to obtain a better description of heavy hadrons, including the $\Lambda$, $\Xi$ and $\Omega$ baryons.

As shown in Fig. 10 the hadronic rescattering has a significant effect on proton $v_{2}$ and we expect a similar effect in the higher harmonics as well. Figures 12 and 13 show our calculations for the identified hadron $v_{3}\{2\}$ and $v_{4}\{2\}$ in several centrality classes. Overall, we find that the level of agreement with data and the effect of rescattering is similar for $v_{3}$ and $v_{4}$ to what was observed for $v_{2}$.

We now turn our attention to Au-Au collisions with $\sqrt{s_{\scriptsize\textrm{NN}}}=200\,\textrm{GeV}$ at RHIC. Figure 14 shows the $p_{T}$-differential $v_{2}$, $v_{3}$, and $v_{4}$ of charged hadrons. The $p_{T}$ spectra and differential $v_{2}$ of identified hadrons are shown in Figs. 15 and 16, respectively. The hybrid approach provides a good description of the charged hadron $v_{n}$ at $p_{T}\lesssim 1$ GeV. At higher $p_{T}$, $v_{2}\{4\}$ and $v_{3}\{2\}$ measurements are overestimated. This is a similar trend as seen at the LHC in Fig. 8. The agreement observed with data for the identified spectra (Fig. 15) is also comparable with the LHC results: calculations describe well the measurements in central collisions but increasing tension is seen in more peripheral bins. As for the identified hadron $v_{2}$ (Fig. 16), we highlight that agreement is distinctly better at RHIC than at the LHC, especially for pions.

A larger tension with experimental data is found for the $p_{T}$ spectra and $p_{T}$-differential elliptic flow of strange baryons, as shown in Fig. 17. We repeat that these heavier hadrons are more sensitive to the transition between hydrodynamics and UrQMD than lighter ones, and that this level of agreement with measurements is not unexpected.

IV Summary and conclusion

In this paper, we compared a hybrid model of IP-Glasma initial conditions, shear and bulk viscous hydrodynamics (MUSIC), and microscopic hadronic transport (UrQMD) with a wide range of integrated and differential measurements from Pb-Pb collisions ($\sqrt{s_{NN}}=2.76\,\textrm{TeV}$) at the LHC and Au-Au collisions ($\sqrt{s_{NN}}=200\,\textrm{GeV}$) at RHIC. We investigated how different observables depend on our model parameters, such as the transport coefficients and the switching temperature from hydrodynamics to the hadronic transport. We found that the bulk viscosity is important to consistently describe the mid-rapidity multiplicity, mean $p_{T}$ of identified hadrons, and the integrated $v_{n}$ within this model.

The inclusion of the bulk viscosity reduces our estimate of the value of the effective shear viscosity by approximately $50\%$. This reduction of shear viscosity is consistent with the intuition that both the shear and bulk viscosities act to reduce the anisotropic flow, and that to produce a similar amount of entropy generated by the larger shear viscosity alone, the shear viscosity in the presence of non-zero bulk viscosity should be smaller.

Heavy hadrons were found to be particularly sensitive to the switching temperature between hydrodynamics and the afterburner. Future improvements on the matching between hydrodynamics and the hadronic transport will be important in reducing this dependence.

It should be emphasized that all three components of our model, IP-Glasma, viscous hydrodynamics and the hadronic after-burner play important roles. The energy deposition mechanism of the IP-Glasma model helps provide a good description of the correct higher flow harmonics, and the large gradient found in the initial energy density enhance the effect of bulk viscosity. The hadronic afterburner is important to improve the description of identified particle observables.

Looking ahead, the addition of mini-jets and jet energy loss will allow us to extend the investigations presented in this work in the intermediate and high $p_{T}$ regions of the observables. Moreover in addition to the observables described above, it is also useful to study the effects of fluctuations and transport coefficients on the event plane correlations and flow harmonics correlations $r_{n}$, which is the subject of a future publication. A further area of possible improvement is the treatment of the non-equilibrium corrections to the thermal distribution functions. So far, our $\delta f$ is species-independent. Making it species-dependent following the line of arguments in Ref. Molnar:2014fva , for example, is an undertaking we leave for future work.

Acknowledgments

This work was supported in part by the Natural Sciences and Engineering Research Council of Canada. SR acknowledges funding of a Helmholtz Young Investigator Group VH-NG-822 from the Helmholtz Association and GSI. JFP was supported in part by the U.S. D.O.E. Office of Science, under Award No. DE-FG02-88ER40388. BPS was supported under DOE contract No. DE-SC0012704 and acknowledges a DOE Office of Science Early Career Award. CG gratefully acknowledges support from the Canada Council for the Arts through its Killam Research Fellowship program. Computations were performed on the Guillimin supercomputer at McGill University under the auspices of Calcul Quebec and Compute Canada. The operation of Guillimin is funded by the Canada Foundation for Innovation (CFI), the National Science and Engineering Research Council (NSERC), NanoQuebec, and the Fonds Quebecois de Recherche sur la Nature et les Technologies (FQRNT). This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02- 05CH11231.

Appendix A Viscous corrections to the momentum distribution of hadrons

At the end of a hydrodynamic simulation, fluid elements must be converted into hadronic degrees of freedom. This conversion is made possible under the assumption that hydrodynamics and kinetic theory have an overlapping region of validity in the late stage of the collision. This overlap allows for the momentum distribution of hadrons to be related to the energy-momentum tensor of the fluid in such a way that energy and momentum are conserved across this transition.

From a kinetic theory point of view, the energy-momentum tensor $T^{\mu\nu}$ only contains information about the second moment of the momentum distribution function, which constrains predominantly the small momentum region of the distribution. In consequence, the transition from fluid to particles carries some ambiguity, since multiple hadronic momentum distributions similar at lower momentum but different at higher momentum can correspond to the same energy-momentum tensor. On the other hand this uncertainty in the higher momentum region of the distribution should not be a major issue for soft hadronic observables, which are the observables of interest in hydrodynamic simulations of heavy ion collisions.

The matching from fluid to hadrons also depends on the collision kernel describing the microscopic interactions of all species of hadrons, which is not known well. Even simplified description of species dependence of hadronic interactions can become quite challenging to handle (see e.g. Ref. Molnar:2014fva ). In consequence, simpler approximations are generally made regarding the collision kernel describing the microscopic interactions of hadrons in order to relate the energy-momentum tensor to the hadron’s momentum distribution. In this work, the relaxation time approximation and the 14-moments approximation were both used to this effect.

An additional assumption made regarding the dependence of the hadronic momentum distribution on the shear stress tensor $\pi^{\mu\nu}$ and the bulk pressure $\Pi$ is that it can be linearized:

[TABLE]

where we used the common notation that the linearized term depending on the shear stress tensor $\pi^{\mu\nu}$ is referred to as $\delta f_{\scriptsize\textrm{shear}}$ and the one depending on the bulk pressure $\Pi$ is $\delta f_{\scriptsize\textrm{bulk}}$. The functional form of $\mathcal{C}_{\scriptsize\textrm{shear}}(P)$ and $\mathcal{C}_{\scriptsize\textrm{bulk}}(P)$ depend on the collision kernel used to describe hadronic interactions. The explicit form of $\delta f_{\scriptsize\textrm{shear/bulk}}$ is given by Eqs (32) and (33).

Since there is a certain level of uncertainty in the determination of $\delta f_{\scriptsize\textrm{shear}}$ and $\delta f_{\scriptsize\textrm{bulk}}$ from the energy-momentum tensor, it is useful to quantify the dependence of hadronic observables on these two quantities. In this section, this is done for the integrated observables shown in Section III.1 for the LHC. Since the effect of the bulk pressure and the shear stress tensor on integrated hadronic observables are significantly different, they are discussed separately in this Appendix.

A.1 Corrections from shear viscosity

Because of the tensor structure of the shear $\pi^{\mu\nu}$-linearized momentum distribution the multiplicity $dN/dy$ does not depend on $\delta f_{\scriptsize\textrm{shear}}$ (for a boost-invariant system). Experiments often employ cuts in transverse momentum when calculating the multiplicity of hadrons, which will lead to some dependence on $\delta f_{\scriptsize\textrm{shear}}$. Nevertheless, these cuts are very small for the measurements used in this work and we verified numerically that the multiplicity of hadrons is essentially identical with and without $\delta f_{\scriptsize\textrm{shear}}$.

The effect of $\delta f_{\scriptsize\textrm{shear}}$ on the average transverse momentum of thermal (Cooper-Frye) pions, kaons and protons and the $p_{T}$-integrated $v_{n}$ of thermal charged hadrons is shown in Figs. 18(a) and (b) respectively. The effect on $\langle p_{T}\rangle$ is very small, and we verified that it remains small even after hadronic decays are taken into account. The $v_{n}$ of charged hadrons displays a larger dependence on $\delta f_{\scriptsize\textrm{shear}}$, of $2-5\%$ for the $v_{2}$, $5-10\%$ for the $v_{3}$ and $10-20\%$ for the $v_{4}$. We verified that we obtain similar numbers after hadronic decays are included.

A.2 Corrections from bulk viscosity

The $\delta f_{\scriptsize\textrm{bulk}}$ used in this work has an explicit dependence on the mass of hadrons — see Eq.(33). Consequently, it is to be expected that the effect of $\delta f_{\scriptsize\textrm{bulk}}$ on different species of hadron will show a mass dependence.

The effect of $\delta f_{\scriptsize\textrm{bulk}}$ on the pion, kaon and proton $dN/dy$ is shown in Fig. 19. Figure 19(a) is for thermal hadrons, and (b) is after hadronic decays. The effect of $\delta f_{\scriptsize\textrm{bulk}}$ decreases the multiplicity of thermal pions by $\sim 15\%$, very slightly decreases the multiplicity of thermal kaons ($\sim 5\%$), and increases the multiplicity of protons by $\sim 10\%$. There is thus a change in the effect of $\delta f_{\scriptsize\textrm{bulk}}$ on the multiplicity at a mass slightly above the kaon mass. Hadronic decays cancel out the suppression from $\delta f_{\scriptsize\textrm{bulk}}$ on thermal pions against the enhancement from $\delta f_{\scriptsize\textrm{bulk}}$ on heavier hadrons which decay into pions. The result is a negligible effect of $\delta f_{\scriptsize\textrm{bulk}}$ on the final pion multiplicity. A similar effect is seen in the final kaon multiplicity, while the already-enhanced thermal proton multiplicity is further increased after decays ($\sim 20\%$) by the effect of $\delta f_{\scriptsize\textrm{bulk}}$ on heavier hadrons.

The effect of $\delta f_{\scriptsize\textrm{bulk}}$ on the average transverse momentum of pions, kaons and protons also shows a mass dependence, as seen in Fig. 20, with (a) being once again for thermal hadrons and (b) being the result including hadronic decays. The thermal pion $\langle p_{T}\rangle$ is suppressed by $\sim 10\%$, while thermal protons are suppressed by $\sim 5\%$ and kaons are in-between. Since heavier hadrons have a smaller correction from $\delta f_{\scriptsize\textrm{bulk}}$, the inclusion of hadronic decays lessens the effect of $\delta f_{\scriptsize\textrm{bulk}}$ on lighter hadrons, as seen in Fig. 20(b).

Finally, unlike for $\delta f_{\scriptsize\textrm{shear}}$, we found that the $\delta f_{\scriptsize\textrm{bulk}}$ used in this work leaves the $p_{T}$-integrated $v_{n}$ of charged hadrons unchanged, whether hadronic decays are included or not.

We highlight that the multiplicity and average transverse momentum of hadrons, which were largely insensitive to $\delta f_{\scriptsize\textrm{shear}}$, are affected by $\delta f_{\scriptsize\textrm{bulk}}$. We thus find the interesting conclusion that integrated observables that are not sensitive to $\delta f_{\scriptsize\textrm{shear}}$ are sensitive to $\delta f_{\scriptsize\textrm{bulk}}$, and vice versa. While the effect of $\delta f_{\scriptsize\textrm{bulk}}$ is not very large, the results found in Figs. 19 and 20 certainly warrants additional investigations in the future about the effect of $\delta f_{\scriptsize\textrm{bulk}}$ on hadronic observables.