3+1d quasiparticle anisotropic hydrodynamics for ultrarelativistic heavy-ion collisions

TL;DR

This paper introduces a 3+1d quasiparticle anisotropic hydrodynamics framework that effectively models the evolution of quark-gluon plasma in heavy-ion collisions, successfully matching experimental data from LHC Pb+Pb collisions.

Contribution

It presents the first application of 3+1d aHydroQP to experimental heavy-ion collision data, incorporating shear and bulk viscosity with a self-consistent freezeout process.

Findings

Good agreement with experimental observables

Effective inclusion of viscous effects

Successful modeling of particle spectra and flow

Abstract

We present the first comparisons of experimental data with phenomenological results from 3+1d quasiparticle anisotropic hydrodynamics (aHydroQP). We compare charged-hadron multiplicity, identified-particle spectra, identified-particle average transverse momentum, charged-particle elliptic flow, and identified-particle elliptic flow produced in LHC 2.76 TeV Pb+Pb collisions. The dynamical equations used for the hydrodynamic stage utilize non-conformal aHydroQP. The resulting aHydroQP framework naturally includes both shear and bulk viscous effects in addition to higher-order non-linear transport coefficients. The 3+1d aHydroQP evolution obtained is self-consistently converted to hadrons using anisotropic Cooper-Frye freezeout performed on a fixed-energy-density hypersurface. The final production and decays of the primordial hadrons are modeled using a customized version of THERMINATOR 2. In this first study, we utilized smooth Glauber-type initial conditions and a single effective freeze-out temperature T_FO = 130 MeV with all hadronic species in full chemical equilibrium. With this rather simple setup, we find a very good description of many heavy-ion observables.

Figures

Equations

$$f(x,p)=f_{\mathrm{iso}}\left(\frac{1}{\lambda }\sqrt{p_{\mu }{\Xi }^{\mu \nu }{p}_{\nu }}\right),$$

$$p^{\mu }{\partial }_{\mu }f+\frac{1}{2}{\partial }_i{m}^2{\partial }_{(p)}^{i}f=-\mathcal{C}[f],$$

Full text

Quark-gluon plasma, Relativistic heavy-ion collisions, Anisotropic hydrodynamics, Equation of state, Boltzmann equation

Ultrarelativistic heavy-ion collision experiments at the Relativistic Heavy Ion Collider and Large Hadron Collider (LHC) were designed to create and study the quark-gluon plasma (QGP). Relativistic hydrodynamics has been quite successful in describing the collective behavior observed in high-energy heavy-ion collisions Huovinen et al. (2001); Hirano and Tsuda (2002); Kolb and Heinz (2003) and the current focus of the relativistic hydrodynamics community is on further improvements of the models to include e.g. bulk viscous effects and higher-order transport coefficients Muronga (2002, 2004); Muronga and Rischke (2004); Heinz et al. (2006); Baier et al. (2006); Romatschke and Romatschke (2007); Song and Heinz (2008a); Dusling and Teaney (2008); Song and Heinz (2008b); Baier et al. (2008); Luzum and Romatschke (2008); Song and Heinz (2009); Heinz (2010); Denicol et al. (2010a, b); Schenke et al. (2011a, b); Bozek (2011); Niemi et al. (2011, 2012); Bożek and Wyskiel-Piekarska (2012); Denicol et al. (2012, 2014a); Ryu et al. (2015); Tinti et al. (2017) (see Heinz and Snellings (2013); Gale et al. (2013); Jeon and Heinz (2016) for recent reviews). The goal of the relativistic viscous hydrodynamics program is to constrain key properties of the QGP such as its initial energy density, initial pressure anisotropies, shear viscosity, bulk viscosity, etc. and to also provide the soft-background evolution necessary to compute QGP-modification of hard probes such as jets and heavy quark bound states.

One of the issues faced by practitioners of traditional second-order viscous hydrodynamics approaches is that, at early times after the nuclear impact, the QGP possesses a high degree of momentum-space anisotropy in the fluid local rest frame, ${\cal P}_{T}/{\cal P}_{L}\gg 1$. The magnitude of the resulting momentum-space anisotropy is large at early times after the initial nuclear impact and also near the transverse/longitudinal “edges” of the QGP at all times. In these spacetime regions, traditional viscous hydrodynamics is being pushed to its limits, resulting in potentially negative total pressures and violations of positivity of the one-particle distribution function Strickland (2014).

As a way to address these problems, it was suggested that one should reorganize the expansion of the one-particle distribution function around a leading-order form which possesses intrinsic momentum-space anisotropies but still guarantees positivity Florkowski and Ryblewski (2011); Martinez and Strickland (2010). This method has become known as anisotropic hydrodynamics (aHydro). Since the two original papers Florkowski and Ryblewski (2011); Martinez and Strickland (2010), there has been a great deal of progress in aHydro Martinez et al. (2012); Ryblewski and Florkowski (2012); Bazow et al. (2014); Tinti and Florkowski (2014); Nopoush et al. (2014); Tinti (2016); Bazow et al. (2015); Strickland et al. (2016); Alqahtani et al. (2015); Molnar et al. (2016a, b); Alqahtani et al. (2017) including applications to cold atomic gases near the unitary limit Bluhm and Schaefer (2015, 2016). In parallel, there have been efforts to construct exact solutions to the Boltzmann equation in some simple cases which can be used to test the efficacy of various dissipative hydrodynamics approaches, and it has been shown that aHydro most accurately reproduces all known exact solutions, even in the limit of very large $\eta/s$ and/or initial momentum-space anisotropy Florkowski et al. (2013a, b); Denicol et al. (2014b, c); Nopoush et al. (2015a); Heinz et al. (2016); Molnar et al. (2016b).

A recent focus of research has been on turning aHydro into a practical phenomenological tool with a realistic equation of state (EoS) and self-consistent anisotropic hadronic freeze-out. In this paper, we present the first comparisons of experimental data with phenomenological results obtained using (1) generalized 3+1d aHydro including three momentum-space anisotropy parameters in the underlying distribution function, (2) the quasiparticle aHydro (aHydroQP) method for implementing a realistic EoS Alqahtani et al. (2015, 2017); Alqahtani and Strickland (2016) and (3) anisotropic Cooper-Frye freezeout Nopoush et al. (2015b); Alqahtani et al. (2017) using the same distribution form as was assumed for the dynamical equations. All previous phenomenological applications of aHydro have relied on the approximate conformal factorization of the energy-momentum tensor, see e.g. Martinez et al. (2012); Ryblewski and Florkowski (2012); Nopoush et al. (2016); Strickland (2016), and/or have used isotropic freezeout Ryblewski and Florkowski (2012). For modeling the primordial hadron production and subsequent hadronic decays we use a customized version of THERMINATOR 2 which has been modified to accept ellipsoidally anisotropic distribution functions Chojnacki et al. (2012).

1. Model: In aHydro, the leading-order one-particle distribution function is assumed to be of generalized Romatschke-Strickland form Romatschke and Strickland (2003); Martinez et al. (2012); Nopoush et al. (2014)

[TABLE]

where $\lambda$ is an energy scale which resembles the temperature in the anisotropic distribution, $\Xi^{\mu\nu}\equiv u^{\mu}u^{\nu}+\xi^{\mu\nu}-\Phi\Delta^{\mu\nu}$ is the anisotropy tensor, $\xi^{\mu\nu}$ obeys $u_{\mu}\xi^{\mu\nu}=0$ and ${\xi^{\mu}}_{\mu}=0$, $\Phi$ is the bulk degree of freedom, and $\Delta^{\mu\nu}=g^{\mu\nu}-u^{\mu}u^{\nu}$ is the transverse projector. Since $\xi^{\mu\nu}$ is traceless and orthogonal to $u^{\mu}$, there are five independent components. In this work, we assume that $\xi^{\mu\nu}$ is diagonal, $\xi^{\mu\nu}={\rm diag}(0,{\boldsymbol{\xi}})$, in which case it has only two independent degrees of freedom. Taken together with $\Phi$, this gives three independent degrees of freedom which map to three ellipsoidal anisotropies in momentum space. It is assumed that the function $f_{\rm iso}(x)$ is a thermal distribution function which can be identified as a Fermi-Dirac, Bose-Einstein, or classical Boltzmann distribution. Herein, we use the Boltzmann distribution when computing the various moment integrals entering the aHydroQP equations of motion and we assume that the fluid chemical potentials are zero at all times. When freezing-out to specific hadron types, however, we use the quantum distribution appropriate for each particle type.

In order to obtain the dynamical equations necessary, we consider a system of quasiparticles with energy-density-dependent masses. In this case, the Boltzmann equation is Jeon and Yaffe (1996)

[TABLE]

with ‘$i$’ labeling spatial indices and ${\cal C}[f]$ being the collisional kernel which, herein, we treat in relaxation-time approximation. In order to conserve energy-momentum and maintain thermodynamic consistency in equilibrium, one must introduce an additional degree of freedom to the energy-momentum tensor, $T^{\mu\nu}=T^{\mu\nu}_{\rm kinetic}+Bg^{\mu\nu}$, where, in general, $B$ is a function of all system parameters and $g^{\mu\nu}$ is the metric tensor Gorenstein and Yang (1995); Jeon and Yaffe (1996); Romatschke (2012); Alqahtani et al. (2015, 2017). By taking momentum-moments of the quasiparticle Boltzmann equation (2), one can obtain a system of partial differential equations for the ellipsoidal anisotropy parameters $\boldsymbol{\xi}$, the scale parameter $\lambda$, and the fluid four-velocity $u^{\mu}$. These dynamical equations form the basis of aHydroQP Alqahtani et al. (2015). The three anisotropies encode the effects of both shear and bulk viscous corrections and, as usual, the fluid four-velocity is normalized to unity.

Using the quasiparticle setup, one can extract the bulk viscosity in the near-equilibrium limit. In Refs. Romatschke (2012) and Tinti et al. (2017) one can find expressions for the bulk viscosity to entropy density $\zeta/s$ in Eqs. (4.4) and (45), respectively. When evaluated, both expressions give the same result for a system of quasiparticles with a temperature-dependent mass. The result is plotted as a black solid line in Fig. 1. In addition to this exact analytic expression, which is valid for all values of $m/T$, we plot an often-used small-mass expansion result, $\zeta/s=15\eta/s\,(1/3-c_{s}^{2})^{2}$, for purposes of comparison as a red dashed line. For both curves we assumed that $\eta/s=0.159$.

In the inset of Fig. 1 we plot the extracted value of $m/T$ obtained by fitting to the Wuppertal-Budapest continuum-extrapolated results for the QCD entropy density Alqahtani et al. (2015). As can be seen from the inset, at small temperatures, the value of $m/T$ necessary to fit the lattice data Borsanyi et al. (2010) is not small, invalidating commonly used small-mass approximations. As Fig. 1 demonstrates, the quasiparticle model used herein has a finite bulk viscosity to entropy density ratio which peaks in the vicinity of the phase transition from QGP to a hadronic gas; however, the magnitude of the peak is much smaller than many other phenomenologically used ansätze for $\zeta/s$. For example, in Ref. Ryu et al. (2015) the authors have a peak value of $\zeta/s$ which is approximately $0.3$.

2. Results and Discussions: In this paper we present comparisons of our aHydroQP results with $\sqrt{s_{NN}}$ = 2.76 TeV Pb+Pb collision data available from the ALICE collaboration. For our initial condition we take the system to be isotropic in momentum space with zero transverse flow and Bjorken flow in the longitudinal direction. In the transverse plane, the initial energy density is computed from a linear combination of smooth Glauber wounded-nucleon and binary-collision profiles with a binary mixing factor of $\alpha=0.15$. In the longitudinal direction, we used a “tilted” profile with a central plateau and Gaussian “wings” resulting in a profile function of the form $\rho(\varsigma)\equiv\exp\left[-(\varsigma-\Delta\varsigma)^{2}/(2\sigma_{\varsigma}^{2})\,\Theta(|\varsigma|-\Delta\varsigma)\right]$, with $\varsigma=\operatorname{arctanh}(z/t)$ being spatial rapidity. The parameters entering the longitudinal profile function were fitted to the pseudorapidity distribution of charged hadrons with the results being $\Delta\varsigma=2.3$ and $\sigma_{\varsigma}=1.6$. The first quantity sets the width of the central plateau and the second sets the width of the Gaussian “wings”.

The resulting initial energy density at a given transverse position ${\bf x}_{\perp}$ and spatial rapidity $\varsigma$ was computed using ${\cal E}\propto(1-\alpha)\rho(\varsigma)\left[W_{A}({\bf x}_{\perp})g(\varsigma)+W_{B}({\bf x}_{\perp})g(-\varsigma)\right]+\alpha\rho(\varsigma)C({\bf x}_{\perp})$, where $W_{A,B}({\bf x}_{\perp})$ is the wounded nucleon density for nuclei $A$ and $B$, $C({\bf x}_{\perp})$ is the binary collision density, and $g(\varsigma)$ is the “tilt function”. The tilt function $g(\varsigma)=0$ if $\varsigma<-y_{N}$, $g(\varsigma)=(\varsigma+y_{N})/(2y_{N})$ if $-y_{N}\leq\varsigma\leq y_{N}$, and $g(\varsigma)=1$ if $\varsigma>y_{N}$ where $y_{N}=\log(2\sqrt{s_{NN}}/(m_{p}+m_{n}))$ is the nucleon momentum rapidity Bozek and Wyskiel (2010).

For all results presented herein, we solved the aHydroQP dynamical equations on a $64^{3}$ lattice with lattice spacings $\Delta x=\Delta y=0.5$ fm and $\Delta\varsigma$ = 0.375. We computed spatial derivatives using fourth-order centered-differences and, for temporal updates, we used fourth-order Runge-Kutta with step size of $\Delta\tau=0.02$ fm/c and a weighted-LAX smoother to regulate potential numerical instabilities associated with the centered-differences scheme Martinez et al. (2012). We started the aHydroQP evolution at $\tau_{0}=0.25$ fm/c and ended it when the highest effective temperature in the three-volume was sufficiently below the freeze-out temperature.

After running the full 3+1d evolution of the system using aHydroQP, we extracted a fixed energy-density freeze-out hypersurface corresponding to a given effective temperature. The fluid anisotropy tensor and scale parameter were assumed to be the same for all hadronic species and we additionally assumed that all produced hadrons were in chemical equilibrium. The aHydroQP distribution function parameters on the freezeout hypersurface were fed into a customized version of THERMINATOR 2 which uses Monte-Carlo sampling to generate final hadronic configurations. Once the primordial hadrons were sampled in this manner, the subsequent hadronic decays proceeded as usual. In the plots shown herein, we used between 7,400 and 36,200 hadronic events, depending on the centrality class and the target observable, e.g. for centrality classes in which we show identified-particle $v_{2}(p_{T})$, more hadronic events were used in order to increase statistics. In all plots, the statistical uncertainty of our model results associated with hadronic Monte-Carlo sampling is indicated by a shaded band surrounding the central line, which indicates the hadronic event-averaged value.

To fix the remaining model parameters, we used scans in the initial central temperature $T_{0}$, the freezeout temperature $T_{\rm FO}$, and $\eta/s$, where the latter was assumed to be a temperature-independent constant. The theoretical predictions resulting from this scan were compared to experimental data from the ALICE collaboration for the differential spectra of pions, kaons, and protons in both the 0-5% and 30-40% centrality classes. The fitting error was minimized across species, with equal weighting for the three particle types. The parameters obtained from this procedure were $T_{0}=600$ MeV, $\eta/s=0.159$, and $T_{\rm FO}=130$ MeV. Herein, $T_{0}$ is the initial temperature which would be obtained in a perfectly central collision at ${\bf x}=0$. The resulting fit to the spectra that emerged in the 0-5% and 20-30% centrality classes is shown in Fig. 2. As can be seen from Fig. 2, the resulting spectra fits are quite good, allowing for a simultaneous description of the pion, kaon, and proton spectra. Note that, for pions, the model slightly underpredicts the pion spectrum at low transverse momentum. This discrepancy is similar to what is observed in other hydrodynamic models.

Using the parameters determined by the procedure outlined above, we then proceeded to calculate other observables. In Fig. 3, we present three panels which show (a) the charged-hadron multiplicity in different centrality classes as a function of pseudorapidity; (b) the average transverse momentum of pions, kaons, and protons as a function of centrality; and (c) the integrated $v_{2}$ for charged hadrons as a function of centrality. In each panel, we compare to data reported by the ALICE collaboration. As can be seen from panel (a), our model is able to describe the charged hadron multiplicity as a function pseudorapidity quite well in all centrality classes. From panel (b) we see that the model is also able to reproduce the average $p_{T}$ of pions, kaons, and protons quite well. The fit quality achieved is similar to Ref. Ryu et al. (2015), however, we note that in our model the peak value of $\zeta/s$ is substantially smaller than what was assumed in Ref. Ryu et al. (2015). Turning to panel (c) we compare our model predictions computed using the geometrical $v_{2}\sim\langle\cos(2\phi)\rangle$ for all charged hadrons with ALICE results obtained using second- and fourth-order cumulants $v_{2}\{2\}$ and $v_{2}\{4\}$. As we can see from panel (c), the model underestimates $v_{2}$ for very central collisions. This is to be expected since we did not include event-by-event fluctuations in the initial condition. That being said, we see that the qualitative behavior of $v_{2}$ as a function of centrality is well reproduced.

Finally, turning to Fig. 4 we present two panels which compare our model predictions for the identified-particle $v_{2}(p_{T})$ with experimental data from the ALICE collaboration. The top and bottom panels show the results obtained in the 20-30% and 30-40% centrality classes, respectively. As can be seen from these panels, the model provides a very good description of the identified-particle elliptic flow. In the 20-30% centrality class, the model is in good agreement with the pion, kaon, and proton data out to $p_{T}\sim$ 1.5, 1.5, and 2.5 GeV, respectively. In the 30-40% centrality class, the model is in good agreement with the pion, kaon, and proton data out to $p_{T}\sim$ 1, 1, and 2 GeV, respectively. In order to improve the agreement between theory and data, it would seem that one has to, at the very least, relax the assumption of a temperature independent $\eta/s$. Since high-momentum hadrons are produced significantly at early times after the collision when the energy density is high, it is natural to expect that elliptic flow would be reduced at high $p_{T}$ since the effective shear viscosity would be larger.

3. Conclusions and Outlook: In this paper, we have provided the first phenomenological comparisons of aHydroQP with LHC experimental data, which includes (1) generalized aHydroQP including three momentum-space anisotropy parameters in the underlying distribution function, (2) the quasiparticle method for implementing a realistic EoS and (3) anisotropic Cooper-Frye freezeout using the same distribution form as was assumed for the dynamical equations. For modeling the primordial hadron production and subsequent hadronic decays we used a customized version of THERMINATOR 2 which has been modified to accept ellipsoidally anisotropic distribution functions. For this initial application, we assumed smooth Glauber initial conditions which were a linear combination of wounded-nucleon and binary-collision profiles. We further assumed that the system was initially isotropic in momentum space. With these assumptions, we performed a parameter scan and, through comparisons of the identified-particle spectra emerging from the model and experimental data, we were able to find a best fit with $T_{0}=600$ MeV, $\eta/s=0.159$, and $T_{\rm FO}=130$ MeV. With this small set of parameters we were able to obtain good agreement between the model and experimental data for the identified-particle spectra, the identified-particle average transverse momentum as a function of centrality, the charged-hadron multiplicity as a function of pseudorapidity, the charged-particle $v_{2}$ as a function of centrality, and the identified-particle $v_{2}$ as a function of transverse momentum. We note, in particular, that we were able to obtain a good description of the average transverse momentum of pions, kaons, and protons with a much smaller peak value for the bulk viscosity to entropy density ratio than previous studies (see e.g. Ref. Ryu et al. (2015)). This suggests that there is a fair amount of hydrodynamical model variation in statements about the magnitude of the bulk viscosity in the QGP.

This study provides the first solid evidence that it is possible to apply aHydroQP to obtain a successful phenomenological description of the QGP. Of course, this study is only the first major step. Looking forward, it is necessary to include realistic fluctuating initial conditions, temperature-dependent shear viscosity to entropy density ratio, realistic initial momentum-space anisotropy profiles (tied to the fluctuating initial conditions), more realistic collisional kernels, etc. These are all currently in progress and, based on the results obtained herein using a somewhat simple setup, we are quite optimistic that aHydroQP can be used as a reliable modeling tool in the future.

Looking beyond the phenomenological applications, we emphasize that the aHydroQP formalism used herein represents an important step forward in the self-consistent implementation of both large shear corrections and non-conformal effects in relativistic dissipative hydrodynamics. Recent works have shown that, in the context of second-order viscous hydrodynamics, the self-consistent incorporation of the temperature-dependence of the quasiparticle mass (vHydroQP) results in important modifications to QGP transport coefficients, particularly at low temperatures Tinti et al. (2017). Additionally, Ref. Tinti et al. (2017) demonstrated that prior aHydroQP results published in Ref. Alqahtani et al. (2015) agree very well with the vHydroQP evolution and that both quasiparticle formalisms result in a qualitatively different evolution of the bulk viscous correction compared to existing approaches which ignore the temperature dependence of the quasiparticle mass. The phenomenological results obtained herein suggest that, if non-conformal aspects are more carefully taken into account, one can obtain a very good description of many key heavy-ion observables.