Anisotropic hydrodynamic modeling of 2.76 TeV Pb-Pb collisions

TL;DR

This paper demonstrates that 3+1d quasiparticle anisotropic hydrodynamics accurately models various observables in 2.76 TeV Pb-Pb collisions at the LHC, matching experimental data well.

Contribution

It introduces a quasiparticle anisotropic hydrodynamics model with temperature-dependent mass fitted to lattice QCD, improving the description of heavy-ion collision dynamics.

Findings

Good agreement with experimental particle spectra.

Accurate reproduction of elliptic flow and HBT radii.

Model captures shear and bulk viscosity effects.

Abstract

We compare phenomenological results from 3+1d quasiparticle anisotropic hydrodynamics (aHydroQP) with experimental data collected in LHC 2.76 TeV Pb-Pb collisions. In particular, we present comparisons of particle spectra, average transverse momentum, elliptic flow, and HBT radii. The aHydroQP model relies on the introduction of a single temperature-dependent quasiparticle mass which is fit to lattice QCD data. By taking moments of the resulting Boltzmann equation, we obtain the dynamical equations used in the hydrodynamic stage which include the effects of both shear and bulk viscosities. At freeze-out, we use anisotropic Cooper-Frye freeze-out performed on a fixed-energy-density hypersurface to convert to hadrons. To model the production and decays of the hadrons we use THERMINATOR 2 which is customized to sample from ellipsoidal momentum-space distribution functions. Using smooth Glauber initial conditions, we find very good agreement with many heavy-ion collision observables.

Figures

Equations

$$u^{\mu }$$

$$X^{\mu }$$

$$Y^{\mu }$$

$$Z^{\mu }$$

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

$${\alpha }_i\equiv (1+{\xi }_i+\Phi {)}^{-1/2},$$

$$f(x,p)=f_{\mathrm{eq}}\left(\frac{1}{\lambda }\sqrt{\sum _i\frac{p_i^2}{{\alpha }_i^2}+m^2}\right).$$

$$\frac{\eta }{{\tau }_{\mathrm{eq}}}=\frac{1}{T}{I}_{3,2}({\hat{m}}_{\mathrm{eq}}),$$

$$I_{3,2}(x)$$

$$K_{i,1}(x)$$

$${\tau }_{\mathrm{eq}}(T)=\bar{\eta }\frac{\mathcal{E}+\mathcal{P}}{I_{3,2}({\hat{m}}_{\mathrm{eq}})},$$

$$\frac{\zeta }{{\tau }_{\mathrm{eq}}}=\frac{5}{3T}{I}_{3,2}({\hat{m}}_{\mathrm{eq}})-c_s^2(\mathcal{E}+\mathcal{P})+c_s^{2}m\frac{dm}{dT}{I}_{1,1}({\hat{m}}_{\mathrm{eq}}),$$

$$I_{1,1}(x)=\frac{N_{\rm dof}T^{3}\,x^{3}}{6\pi^{2}}\Bigg{[}\frac{1}{4}\Big{(}K_{3}(x)-5K_{1}(x)\Big{)}+K_{i,1}(x)\Bigg{]}\,.$$

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

$$T^{\mu \nu }=T_{\mathrm{kinetic}}^{\mu \nu }+B(T)g^{\mu \nu },$$

$$T^{\mu \nu }=\mathcal{E}{u}^{\mu }{u}^{\nu }+{\mathcal{P}}_x{X}^{\mu }{X}^{\nu }+{\mathcal{P}}_y{Y}^{\mu }{Y}^{\nu }+{\mathcal{P}}_z{Z}^{\mu }{Z}^{\nu }.$$

$$\mathcal{E}$$

$${\mathcal{P}}_i$$

$$D_u\mathcal{E}+\mathcal{E}{\theta }_u+{\mathcal{P}}_x{u}_{\mu }{D}_x{X}^{\mu }+{\mathcal{P}}_y{u}_{\mu }{D}_y{Y}^{\mu }+{\mathcal{P}}_z{u}_{\mu }{D}_z{Z}^{\mu }$$

$$D_x{\mathcal{P}}_x+{\mathcal{P}}_x{\theta }_x-\mathcal{E}{X}_{\mu }{D}_u{u}^{\mu }-{\mathcal{P}}_y{X}_{\mu }{D}_y{Y}^{\mu }-{\mathcal{P}}_z{X}_{\mu }{D}_z{Z}^{\mu }$$

$$D_y{\mathcal{P}}_y+{\mathcal{P}}_y{\theta }_y-\mathcal{E}{Y}_{\mu }{D}_u{u}^{\mu }-{\mathcal{P}}_x{Y}_{\mu }{D}_x{X}^{\mu }-{\mathcal{P}}_z{Y}_{\mu }{D}_z{Z}^{\mu }$$

$$D_z{\mathcal{P}}_z+{\mathcal{P}}_z{\theta }_z-\mathcal{E}{Z}_{\mu }{D}_u{u}^{\mu }-{\mathcal{P}}_x{Z}_{\mu }{D}_x{X}^{\mu }-{\mathcal{P}}_y{Z}_{\mu }{D}_y{Y}^{\mu }$$

$${\mathcal{I}}^{\mu \nu \lambda }\equiv \int dP{p}^{\mu }{p}^{\nu }{p}^{\lambda }f(x,p).$$

$$\mathcal{I}$$

$${\mathcal{I}}_i$$

$${\mathcal{I}}_{\mathrm{eq}}(\lambda ,m)$$

$$D_u{\mathcal{I}}_x+{\mathcal{I}}_x({\theta }_u+2u_{\mu }{D}_x{X}^{\mu })$$

$$D_u{\mathcal{I}}_y+{\mathcal{I}}_y({\theta }_u+2u_{\mu }{D}_y{Y}^{\mu })$$

$$D_u{\mathcal{I}}_z+{\mathcal{I}}_z({\theta }_u+2u_{\mu }{D}_z{Z}^{\mu })$$

$${\mathcal{H}}_3(\mathbf{\alpha },\hat{m}){\lambda }^4={\mathcal{H}}_{3,\mathrm{eq}}(1,{\hat{m}}_{\mathrm{eq}})T^4.$$

$$\bigg{(}p^{0}\frac{dN}{d^{3}p}\bigg{)}_{i}=\frac{{\cal N}_{i}}{(2\pi)^{3}}\int\!f_{i}(x,p)\,p^{\mu}d^{3}\Sigma_{\mu}\,,$$

$$\rho (\varsigma )\equiv \exp \left[-(\varsigma -\Delta \varsigma {)}^2/(2{\sigma }_{\varsigma }^2)\Theta (|\varsigma |-\Delta \varsigma )\right].$$

$${\cal E}({\bf x}_{\perp},\varsigma)\propto(1-\chi)\rho(\varsigma)\Big{[}W_{A}({\bf x}_{\perp})g(\varsigma)+W_{B}({\bf x}_{\perp})g(-\varsigma)\Big{]}+\chi\rho(\varsigma)C({\bf x}_{\perp})\,,$$

$$\displaystyle g(\varsigma)=\left\{\begin{array}[]{lcccc}0&&\mbox{if}&&\varsigma<-y_{N}\,,\\ (\varsigma+y_{N})/(2y_{N})&&\mbox{if}&&-y_{N}\leq\varsigma\leq y_{N}\,,\\ 1&&\mbox{if}&&\varsigma>y_{N}\,,\end{array}\right.\,\,\,\,\,\,\,\,\,\,\,\,\,$$

$$\mathcal{D}$$

$$\tilde{\mathcal{D}}$$

Full text

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

I Introduction

Relativistic hydrodynamics has been quite successful in describing the soft hadron spectra ($p_{T}\lesssim 2$ GeV) and collective flow observed in ultrarelativistic heavy-ion collision (URHIC) experiments at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) Huovinen et al. (2001); Hirano and Tsuda (2002); Kolb and Heinz (2003); 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); Broniowski et al. (2009); Heinz (2010); Denicol et al. (2010a, b); Schenke et al. (2011a); Bozek et al. (2011); Schenke et al. (2011b); Bozek (2011); Niemi et al. (2011, 2012); Bożek and Wyskiel-Piekarska (2012); Bozek and Broniowski (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); Jaiswal and Roy (2016) for recent reviews). Early on, the success of this program suggested that the quark-gluon plasma (QGP) generated in URHICs was nearly isotropic in the local rest frame (LRF); however, in practice, one finds that there are rather large momentum-space anisotropies, which are driven primarily by the rapid longitudinal expansion of the QGP created in URHICs Ryblewski (2013); Strickland (2014). All studies indicate that at early times after the nuclear impact, the QGP possesses a high degree of momentum-space anisotropy in the fluid LRF, ${\cal P}_{T}/{\cal P}_{L}\gg 1$, which only slowly relaxes towards unity during the QGP liftetime. Additionally, at all proper times, there are large momentum-space anisotropies near the transverse/longitudinal “edges” of the fireball where the system is nearly free streaming. The situation only gets worse as one goes from AA to pA and pp collisions, since gradients are larger and system lifetimes are considerably shorter. As a consequence, larger non-equilibrium deviations are expected for these systems, which push traditional viscous hydrodynamics to its limits Nopoush et al. (2015a). This has motivated the investigation of alternative formulations of dissipative relativistic hydrodynamics which can be applied to systems which might possess a high degree of momentum-space anisotropy at all points in spacetime Ryblewski (2013); Strickland (2014).

One way to proceed is to 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) and there are now many groups pursuing this idea Ryblewski and Florkowski (2011); 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. (2017a); Bluhm and Schaefer (2015, 2016). One of the selling points for the aHydro approach has been that it better reproduces exact solutions to the Boltzmann equation compared to traditional near-equilibrium viscous hydrodynamic approaches, even in the limit of very large shear viscosity to entropy density ratio and/or initial momentum-space anisotropy Florkowski et al. (2013a, b); Denicol et al. (2014b, c); Nopoush et al. (2015b); Heinz et al. (2015); Molnar et al. (2016b); Martinez et al. (2017). Given this success, the focus has recently turned to making aHydro a practical phenomenological tool with a realistic equation of state (EoS) and self-consistent anisotropic hadronic freeze-out. In a previous short paper Alqahtani et al. (2017b), we presented the first comparisons of experimental data with phenomenological results obtained using generalized 3+1d aHydro including: (1) 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, 2017a); Alqahtani and Strickland (2016), and (3) anisotropic Cooper-Frye freeze-out Nopoush et al. (2015a); Alqahtani et al. (2017a) using the same form for the distribution function as was assumed for the dynamical equations.

To the best of our knowledge, all previous phenomenological applications of aHydro have either 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. (2017); Strickland (2016), and/or have used isotropic freeze-out Ryblewski and Florkowski (2012). To address the issue with freeze-out, we use a customized version of THERMINATOR 2 Chojnacki et al. (2012) which has been modified to accept ellipsoidally-anisotropic distribution functions. In this paper, we extend our previous work Alqahtani et al. (2017b) by (1) giving the details of the formalism of 3+1d aHydroQP and (2) presenting comparisons with a wider set of heavy-ion observables. Here we present comparisons of charged-hadron multiplicity, identified-particle spectra, identified-particle average transverse momentum, charged-particle elliptic flow, identified-particle elliptic flow, the integrated elliptic flow vs pseudorapidity, and the HBT radii. For some observables, such as the spectra, compared to Ref. Alqahtani et al. (2017b) we present comparisons in more centrality classes and with higher statistics.

The structure of the paper is as follows. In Sec. II, we present our general setup. In Sec. III, we introduce the formalism for quasiparticle anisotropic hydrodynamics and then derive the 3+1d dynamical equations. In Sec. IV, we discuss anisotropic Cooper-Frye freeze-out. In Sec. V, we compare our model results obtained using the 3+1d aHydroQP model for Pb-Pb collisions at LHC energies with data from the ALICE collaboration. Sec. VI contains our conclusions and an outlook for the future. In App. A, we list the derivatives used in the body of the paper. App. B contains details concerning the thermodynamic integrals introduced in the body of the text and our optimized scheme for their evaluation.

II Setup

II.1 Conventions and notation

Here we define some conventions and notation that will be used in the body of this paper. The metric is taken to be “mostly minus” with $x^{\mu}=(t,x,y,z)$ where the line element is $ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=dt^{2}-dx^{2}-dy^{2}-dz^{2}$ with $g^{\mu\nu}$ being metric tensor in Minkowski space. The longitudinal proper time is $\tau=\sqrt{t^{2}-z^{2}}$ and the longitudinal spacetime rapidity is $\varsigma=\tanh^{-1}(z/t)$.

The basis vectors for a general 3+1d system in the laboratory frame can be obtained by the following parametrization, which is based on a set of Lorentz transformations applied to the LRF basis vectors Martinez et al. (2012); Ryblewski and Florkowski (2011). The set of successive transformations correspond to a longitudinal boost $\vartheta$ along the beam line, a rotation $\varphi\equiv\tan^{-1}(u_{y}/u_{x})$ around the beam line, and a transverse boost $\theta_{\perp}$, which together yield

[TABLE]

where $u^{\mu}$ is the fluid four-velocity and $u_{\perp}\equiv\sqrt{u_{x}^{2}+u_{y}^{2}}=\sqrt{u_{0}^{2}-1}=\sinh\theta_{\perp}$.

II.2 Distribution Function

The leading order distribution function in aHydro is assumed to be of generalized Romatschke-Strickland form Romatschke and Strickland (2003, 2004)

[TABLE]

where $\lambda$ is an energy scale which becomes the temperature in the isotropic equilibrium limit. The anisotropy tensor has the form $\Xi^{\mu\nu}\equiv u^{\mu}u^{\nu}+\xi^{\mu\nu}-\Delta^{\mu\nu}\Phi$ where $\xi^{\mu\nu}$ is a symmetric traceless tensor obeying $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 Martinez et al. (2012); Nopoush et al. (2014). Using the ellipsoidal form (2) and the tracelessness of $\xi^{\mu\nu}$, we are left with three independent parameters out of the four original parameters $\Phi$ and ${\boldsymbol{\xi}}=(\xi_{x},\xi_{y},\xi_{z})$. In thermal equilibrium, the distribution function $f_{\rm eq}(x)$ can be identified as Fermi-Dirac, Bose-Einstein, or Maxwellian distribution. Herein, ignoring the quantum statistics, we take the Boltzmann form with zero chemical potential.

We note that, in order to perform the integrals, it is useful to define three parameters $\alpha_{i}$ as the system’s independent anisotropy parameters

[TABLE]

such that the distribution function can be written as

[TABLE]

II.3 The equation of state

We use an analytic parameterization of lattice QCD (LQCD) data for the trace anomaly taken from the Wuppertal-Budapest collaboration Borsanyi et al. (2010). This analytic parameterization can be used to calculate the energy density, pressure, and entropy density using standard thermodynamic identities. In order to implement the equation of state in the quasiparticle model, we fit a single temperature-dependent mass, $m(T)$, to the LQCD entropy density. For details, we refer the reader to Ref. Alqahtani et al. (2015).

II.4 Shear and bulk viscosity in the quasiparticle (QP) model

The relaxation time $\tau_{\rm eq}$ can be related to the shear viscosity, $\eta$, for the system of quasiparticles. The shear viscosity for a quasiparticle gas can be found in Refs. Romatschke (2012) and Tinti et al. (2017), Eqs. (4.3) and (46), respectively, with both giving the same result

[TABLE]

with $\hat{m}_{\rm eq}\equiv m/T$ and

[TABLE]

where $N_{\rm dof}$ is the number of degrees of freedom, $K_{n}$ are modified Bessel functions of the second kind, and ${\cal S}_{n}$ are modified Struve functions. Using Eq. (5), the relaxation time can be written as

[TABLE]

where $\bar{\eta}\equiv\eta/s$ with $s$ being the entropy density, ${\cal E}$ is the energy density, and ${\cal P}$ is the pressure.

Using the quasiparticle model, one can extract the bulk viscosity in a similar manner. Expressions for the bulk viscosity for a quasiparticle gas can be found in Refs. Romatschke (2012) and Tinti et al. (2017), Eqs. (4.4) and (45), respectively, with, both again giving the same result

[TABLE]

where

[TABLE]

In Fig. 1-a, we plot the result for the bulk viscosity to entropy density ratio $\zeta/s$ in the quasiparticle model as a black solid line. For comparison, we plot $\zeta/s=15\eta/s\,(1/3-c_{s}^{2})^{2}$, where $c_{s}$ is the speed of sound. This is an often-used small-mass expansion result Weinberg (1971). For both curves, we assume that $\eta/s=2/(4\pi)$.

In Fig. 1-b we plot $m/T$ obtained by fitting to the Wuppertal-Budapest LQCD results for the entropy density Alqahtani et al. (2015). We point out that, at small temperatures, the value of $m/T$ necessary to fit the LQCD data Borsanyi et al. (2010) is not small, invalidating the often used small $m/T$ expansion used to compute effective transport coefficients. The quasiparticle model has a finite bulk viscosity to entropy density ratio as shown in Fig. 1-a. It peaks in the vicinity of the phase transition temperature from QGP to a hadronic gas. By comparing to other ansätze for $\zeta/s$ used in other studies we see that, in our quasiparticle model, the peak value is much smaller $\sim 0.05$, compared to prior works. For example, in Ref. Ryu et al. (2015) the respective peak value of $\zeta/s$ is approximately $0.3$.

III dynamical equations

For a system of quasiparticles with a temperature-dependent mass, the Boltzmann equation is Jeon and Yaffe (1996); Romatschke (2012); Alqahtani et al. (2015)

[TABLE]

with $i\in\{x,y,z\}$. Herein, we take the collisional kernel ${\cal C}[f]$ in the relaxation time approximation, $C[f]=p_{\mu}u^{\mu}(f-f_{\rm eq})/\tau_{\rm eq}$. By taking moments of this equation we can generate the necessary dynamical equations.

The kinetic part of the energy-momentum tensor can be obtained from the second moment of distribution function, however, when the quasiparticle mass is temperature dependent, this quantity is not conserved by itself. In order to enforce thermodynamic consistency and energy-momentum conservation, one must introduce the background field contribution Gorenstein and Yang (1995), $B(T)$, which is fixed through comparison with LQCD data Romatschke (2012); Alqahtani et al. (2015)

[TABLE]

where $T^{\mu\nu}_{\rm kinetic}=\int dP\,p^{\mu}p^{\nu}f(x,p)$ with $dP\equiv E^{-1}\,d^{3}p/(2\pi)^{3}$ being the Lorentz-invariant integration measure.

For the case considered here, namely a diagonal anisotropy tensor, the full energy-momentum tensor can be expanded as

[TABLE]

The resulting energy density and pressures are

[TABLE]

with $i\in\{x,y,z\}$ and the ${\cal H}$-functions appearing above defined in App. B.

Taking the first moment of Boltzmann equation we obtain four equations

[TABLE]

The second moment of Boltzmann equation involves a rank-3 tensor ${\cal I}$ which is the third moment of the distribution function

[TABLE]

Expanding ${\cal I}$ over the basis vectors one has

[TABLE]

with Nopoush et al. (2014)

[TABLE]

where $\alpha=\prod_{i}\alpha_{i}$. Since we have eight independent variables, $\alpha_{x}$, $\alpha_{y}$, $\alpha_{z}$, $u_{x}$, $u_{x}$, $\vartheta$, $T$, and $\lambda$, we need eight equations to solve the full 3+1d aHydroQP system. With ten equations obtained from the second moment of Boltzmann equation, the system is overdetermined. Therefore, we have to come up with a selection rule to choose a subset of these equations in order to close the set of dynamical equations. The final equations are taken from the three diagonal projections of the equation of motion of the third moment, $X_{\mu}X_{\nu}\partial_{\alpha}{\cal I}^{\alpha\mu\nu}$, $Y_{\mu}Y_{\nu}\partial_{\alpha}{\cal I}^{\alpha\mu\nu}$, and $Z_{\mu}Z_{\nu}\partial_{\alpha}{\cal I}^{\alpha\mu\nu}$ giving Alqahtani et al. (2015)

[TABLE]

Finally, in order to compute the local effective temperature $T$, we match the non-equilibrium kinetic energy density with the equilibrium kinetic energy density

[TABLE]

To summarize, to perform the numerical simulations reported herein we use the eight equations resulting from the first (14) and second (18) moments of the Boltzmann equation, together with the matching condition (19).

IV Anisotropic freeze-out

The QGP undergoes freeze-out at late times/low temperatures and the degrees of freedom need to be changed from hydrodynamical variables to hadronic positions and momenta. In this work, we perform “anisotropic Cooper-Frye freeze-out” using Eq. (4) as the form for the one-particle distribution function. The anisotropic distribution function used in the freeze-out is guaranteed to be positive-definite, by construction, in all regions in phase space, avoiding the usual problems encountered within standard viscous hydrodynamic freeze-out. In practice, we construct a constant energy-density hypersurface, defined through $T_{\rm FO}={\cal E}^{-1}({\cal E}_{\rm FO})$. Then, by computing the number of particles that cross this hypersurface, one can determine the number of hadrons produced in heavy-ion collisions at freeze-out using

[TABLE]

where $i$ labels the hadronic species, ${\cal N}_{i}\equiv 2s_{i}+1$ is the degeneracy factor with $s_{i}$ being the spin of particle species $i$, and $f_{i}$ is the distribution function for particle species $i$ taking into account the appropriate quantum statistics.111In THERMINATOR 2, different isospin states are treated separately negating the need for an explicit isospin degeneracy factor. For more details we refer the reader to Nopoush et al. (2015a); Alqahtani et al. (2017a).

V Numerical results

In this section, we present comparisons of our aHydroQP model results with $\sqrt{s_{NN}}$ = 2.76 TeV Pb-Pb collision data available from the ALICE collaboration. To set the initial conditions, we assume the system to be initially isotropic in momentum space ($\alpha_{i}(\tau_{0})=1$), with zero transverse flow (${\bf u}_{\perp}(\tau_{0})=0$), and Bjorken flow in the longitudinal direction ($\vartheta(\tau_{0})=\eta$). The initial energy density distribution in the transverse plane is computed from a “tilted” profile Bozek and Wyskiel (2010). The distribution used is a linear combination of smooth Glauber wounded-nucleon and binary-collision density profiles, with a binary-collision mixing factor of $\chi=0.15$. In the longitudinal direction, we used a profile with a central plateau and Gaussian “tails”, resulting in a longitudinal profile function of the form

[TABLE]

The parameters entering (21) 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 “tails”.

The resulting initial energy density at a given transverse position ${\bf x}_{\perp}$ and spatial rapidity $\varsigma$ was computed using

[TABLE]

where $W_{A,B}({\bf x}_{\perp})$ is the wounded nucleon density for nucleus $A$ or $B$, $C({\bf x}_{\perp})$ is the binary collision density, and $g(\varsigma)$ is the “tilt function”. The tilt function is defined through

[TABLE]

where $y_{N}=\log(2\sqrt{s_{NN}}/(m_{p}+m_{n}))$ is the nucleon momentum rapidity Bozek and Wyskiel (2010).

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. To compute spatial derivatives we used fourth-order centered-differences and, for temporal updates, we used fourth-order Runge-Kutta with step size of $\Delta\tau=0.02$ fm/c. To regulate potential numerical instabilities associated with the centered-differences scheme, we used a weighted-LAX smoother Martinez et al. (2012). In most cases, we set the weighted-LAX fraction to be $0.005$, however, for large impact parameters we used $0.02$.222This does not affect the evolution considerably since, for high impact parameters, the system reaches $T_{\rm FO}$ at times $\lesssim 4$ fm/c. The aHydroQP evolution was started at $\tau_{0}=0.25$ fm/c and stopped when the highest effective temperature in the entire volume was sufficiently below $T_{\rm FO}$.

Using aHydroQP, we first ran the full 3+1d evolution of the system, then we extracted a freeze-out hypersurface based on the effective temperature. We assumed that all hadronic species were in chemical equilibrium and had the same fluid anisotropy tensor ($\Xi_{\mu\nu}$) and scale parameter ($\lambda$). The distribution function parameters on the freeze-out hypersurface were fed into a customized version of THERMINATOR 2 which allows for an ellipsoidal distribution function of the form given in Eq. (4). THERMINATOR 2 performs sampled event-by-event hadronic production from the exported freeze-out hypersurface using Monte-Carlo sampling. It then performs hadronic feed down (resonance decays) for each sampled event. Depending on the observables under consideration and the centrality class considered, one may need to generate more hadronic events for the purposes of improved statistics. For all plots shown herein we used between 7,400 and 36,200 hadronic events per centrality class. We indicate the statistical uncertainty of our model results associated with the hadronic Monte-Carlo sampling by a shaded band surrounding the hadronic event-averaged value (the central line).

In our model we have three remaining free parameters: (1) the initial central temperature $T_{0}$ obtained in a perfectly central collision at ${\bf x}_{\perp}=0$ and $\varsigma=0$, (2) the freeze-out temperature $T_{\rm FO}$, and (3) $\eta/s$ which is assumed to be a (temperature-independent) constant. In order to fix these parameters we scanned over them and compared the theoretical predictions resulting from this scan with 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 are $T_{0}=600$ MeV, $\eta/s=0.159$, and $T_{\rm FO}=130$ MeV.

We first present our comparisons of the transverse momentum spectra of $\pi^{\pm}$, $K^{\pm}$, and $p+\bar{p}$ in six centrality classes 0-5$\%$, 5-10$\%$, 10-20$\%$, 20-30$\%$, 30-40$\%$, and 40-50$\%$ in Fig. 2. These comparisons show that our model provides a very good simultaneous description of the ALICE data for the pion, kaon, and proton spectra Abelev et al. (2013), with largest differences at $p_{T}\gtrsim 1.5$ GeV and relatively high centrality classes 30-40$\%$, and 40-50$\%$. We note that our model slightly underpredicts the pion spectrum at low transverse momentum which is similar to what is observed in other hydrodynamic models (see e.g. Ref. Ryu et al. (2015)). One possible explanation for this discrepancy that has been suggested is pion condensation Begun et al. (2014).

In Fig. 3, we show the charged-hadron multiplicity in different centrality classes as a function of pseudorapidity, $\eta$. In panel (a), we show the 0-5$\%$, 5-10$\%$, 10-20$\%$, 20-30$\%$, and 30-40$\%$ centrality classes, and in panel (b) we show the 40-50$\%$, 50-60$\%$, 70-80$\%$, 80-90$\%$, and 90-100$\%$ centrality classes. As can be seen from both panels, our model is able to describe the charged hadron multiplicity as a function of pseudorapidity Abbas et al. (2013); Adam et al. (2016a) quite well in all centrality classes. Another observable to consider is the average transverse momentum of pions, kaons, and protons as a function of centrality. This is shown in Fig. 4, where our model is again able to reproduce the data reasonably well.

Next, in Fig. 5, we show the integrated elliptic flow coefficient $v_{2}$ for charged hadrons as a function of centrality. Our model predictions were computed using the geometrical definition of the elliptic flow coefficient, $v_{2}\sim\langle\cos(2\phi)\rangle$, for all charged hadrons. The experimental data were obtained using second- and fourth-order cumulants $v_{2}\{2\}$ and $v_{2}\{4\}$ Abelev et al. (2014). From this figure we see that our model agrees well with $v_{2}\{4\}$ measurements at low centrality, but agrees better with $v_{2}\{2\}$ at higher centrality. One would expect better agreement with $v_{2}\{4\}$ than $v_{2}\{2\}$, since the former has non-flow effects subtracted. The fact that we agree better with $v_{2}\{2\}$ at high centrality could be due to the fact that our smooth initial condition is too simple or that we have not included the off-diagonal components of the anisotropy tensor in the evolution and freeze-out.

In Fig. 6, we present comparisons of the identified-particle $v_{2}$ as a function of $p_{T}$ obtained using our model with experimental data reported by ALICE collaboration Abelev et al. (2015). Our model provides a quite reasonable description of the identified-particle elliptic flow as can be seen in panels (b) and (c), 20-30% and 30-40% centrality classes, respectively. In panel (b), the 20-30% centrality class, we see that our model reproduces the data very well for the pion, kaon, and proton data out to $p_{T}\sim$ 1.5, 1.5, and 2.5 GeV, respectively. A very similar agreement is seen in panel (c), the 30-40% centrality class, where the model is in good agreement with the pion, kaon, and proton data out to $p_{T}\sim$ 1, 1, and 2 GeV, respectively. However, in panels (a) and (d), 10-20% and 40-50% centrality classes, respectively, we see less agreement than panels (b) and (c). For example, we underpredict the pion elliptic flow in the 10-20$\%$ centrality class as can be seen from panel (a). Again, as in the case of Fig. 5, this is related to our use of smooth Glauber initial conditions.

In order to further examine how well our model describes various observables, we look at the pseudorapidity dependence of $v_{2}$ for different centrality classes in Fig. 7. As can be seen from Fig. 7 our model results do not fall fast enough at large pseudorapidity compared to the experimental data Adam et al. (2016b). One possible remedy for this may be including temperature-dependent $\eta/s$, since this has been shown to improve agreement with this observable in the context of viscous hydrodynamics Denicol et al. (2016).

In Fig. 8 we compare the HBT radii predicted by aHydroQP with experimental data from the ALICE collaboration Graczykowski (2014). To compute the HBT radii, we used exactly the same parameters used to describe other observables including the same number of hadronic events in each centrality class. In this set of figures, in the left, middle, and right panels we show $R_{\rm out}$, $R_{\rm side}$, and $R_{\rm long}$, respectively, as a function of the mean transverse momentum of the pair $\pi^{+}\pi^{+}$ in four different centrality classes, 0-5%, 5-10%, 10-20%, and 20-30%. From the left column of this set of figures, we see that our model reproduces the data quite well for $R_{\rm out}$ out to $k_{T}\sim$ 0.6 GeV. In the middle column, we present comparisons of $R_{\rm side}$ where our model shows a good agreement out to $k_{T}\sim$ 0.9 GeV. Lastly, in the right column we compare results for $R_{\rm long}$ which show poorer agreement with the data when compared to the first two columns. However, in most cases, our model predictions are within the error bars of the experimental data, with the biggest differences at low $k_{T}\sim$ 0.2 GeV. This is opposite to what we see in $R_{\rm out}$ and $R_{\rm side}$ where we observe good agreement with the experimental data at low $k_{T}$.

For more comparisons, we compare also the ratios of the HBT radii in Fig. 9. In this set of figures, in the left, middle, and right panels we show $R_{\rm out}/R_{\rm side}$, $R_{\rm out}/R_{\rm long}$, and $R_{\rm side}/R_{\rm long}$, respectively, as a function of the mean transverse momentum of the pair $\pi^{+}\pi^{+}$ in four different centrality classes, 0-5%, 5-10%, 10-20%, and 20-30%. We see that our model was able to reproduce the data quite well in all three panels for all centrality classes shown here.

VI Conclusions and Outlook

In this paper we presented phenomenological comparisons of aHydroQP with LHC experimental data collected in 2.76 TeV Pb-Pb collisions. This work is an extension of a previous letter Alqahtani et al. (2017b). Herein, we gave more details about the formalism used and presented a more thorough comparison between our model and LHC data for a variety of observables. In aHydroQP, we included three momentum anisotropy parameters in the underlying distribution function, both in the dissipative hydrodynamic stage and at freeze-out. We also used a quasiparticle implementation of the LQCD EoS in order to take into account the non-conformality of the system. At freeze-out, we used a customized version of THERMINATOR 2 which was modified to accept anisotropic distribution functions of generalized Romatschke-Strickland form. As a first test, in this work, we used smooth Glauber initial conditions which were obtained from a linear combination of wounded-nucleon and binary-collision profiles. We additionally assumed the system to be initially isotropic in momentum space with no initial transverse flow.

To fix the remaining phenomenological parameters, we performed a parameter scan and compared our results with experimentally observed identified-particle spectra in the 5-10% and 30-40% centrality classes. The resulting set of best fit parameters was $T_{0}=600$ MeV, $\eta/s=0.159$, and $T_{\rm FO}=130$ MeV. After this fitting was complete, we computed an array of different heavy-ion observables, finding quite good agreement between our model and experimental data despite our simple smooth initial condition. We looked at particle multiplicity and spectra, average transverse momentum, $v_{2}$, and HBT radii. Compared to Ref. Alqahtani et al. (2017b), we have added additional centrality classes in some cases and increased the statistics associated with the hadronic Monte-Carlo sampling where necessary.

Combined with what was reported in Alqahtani et al. (2017b), the phenomenological results presented herein represent the first aHydro results to include three separate anisotropy parameters together with the quasiparticle method for imposing the EoS and self-consistent anisotropic freeze-out. Compared to prior results which used a single anisotropy parameter and/or an approximate conformal-factorization implementation of the equation of state Ryblewski and Florkowski (2012); Nopoush et al. (2017); Strickland (2016) we see much better agreement with the pion, kaon, and proton spectra and, relatedly, the total multiplicity as a function of pseudorapidity. Prior studies which used the approximate conformal-factorization implementation of the equation of state dramatically underestimated the low $p_{T}$ spectra Nopoush et al. (2017); Strickland (2016), making this the first phenomenological study within the context of aHydro which is able to reproduce both the experimentally observed spectra and elliptic flow.

Looking to future, there is certainly room for improvements in our model. For example, we are working on including a temperature-dependent shear viscosity to entropy density ratio since in this study it was assumed to be constant. Based on prior studies in the context of viscous hydrodynamics Denicol et al. (2016), there is some hope that this will improve the agreement between our model and the experimental data, in particular, with regards to the pseudorapidity dependence of $v_{2}$. We also plan to include realistic fluctuating initial conditions, realistic initial momentum anisotropy profiles, and more realistic collisional kernels. Additionally, we are planning to look at different collision energies, e.g. RHIC 200 GeV collisions and LHC 5.023 TeV collisions, and different colliding systems, e.g. pA and pp, in the near future. The application of aHydro to pA and pp is of particular interest, since in these systems viscous hydrodynamics is being pushed to its limits, especially at freeze-out Nopoush et al. (2015a).

Finally, it would be interesting to apply this formalism also to even lower energy collisions where it is critical to take into account the finite net baryon density, heat flow, etc. In this context, it would also be interesting to extend the formalism to multicomponent fluids, e.g. two- and three-fluid models, similar to what has been done by the Los Alamos Amsden et al. (1978); Clare and Strottman (1986), Kurchatov Institute Mishustin et al. (1988, 1991); Ivanov et al. (2002), Frankfurt Katscher et al. (1993); Dumitru et al. (1995); Brachmann et al. (1997); Reiter et al. (1998); Brachmann et al. (2000), and GSI Russkikh et al. (1994); Ivanov et al. (1991); Ivanov and Russkikh (1994, 1995); Russkikh and Ivanov (1995); Ayik et al. (1994) groups. Along these lines, there have already been prior anisotropic hydrodynamics studies which have considered multi-component fluids with explicit quark and gluon components Florkowski et al. (2013c, 2015); however, to the best of our knowledge there have thus far not been any attempts to do this in the context of multicomponent hadronic fluids.

Acknowledgements.

We thank Piotr Bożek for useful comments. M. Alqahtani was supported by a PhD fellowship from the Imam Abdulrahman Bin Faisal University, Saudi Arabia. M. Nopoush and M. Strickland were supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award No. DE-SC0013470. R. Ryblewski was supported by the Polish National Science Center grant No. DEC-2012/07/D/ST2/02125.

Appendix A Explicit formulas for derivatives

In this section, we introduce the notation used in our formulation of the general moment-based hydrodynamics equations. Using the definitions

[TABLE]

and four-vectors defined in Eq. (1) one obtains

[TABLE]

The divergences can be defined as

[TABLE]

where $\varphi=\tan^{-1}(u_{y}/u_{x})$. Finally, we list the non-vanishing contractions appearing in the second moment equations

[TABLE]

where $\alpha\in\{u,x,y,z\}$. Note that, using the orthogonality of the basis vectors, i.e. $D_{\alpha}(X^{\mu}u_{\mu})=0$, contractions such as $X^{\mu}D_{\alpha}u_{\mu}$ can be related to the contractions listed above, e.g. $X^{\mu}D_{\alpha}u_{\mu}=-u_{\mu}D_{\alpha}X^{\mu}$.

Appendix B Special functions and derivatives

The ${\cal H}$-functions appearing in the body of the paper can be written as

[TABLE]

where $\boldsymbol{\alpha}=(\alpha_{x},\alpha_{y},\alpha_{z})$, $\hat{m}=m/\lambda$, $i\in\{x,y,z\}$, $\alpha\equiv\prod_{i}\alpha_{i}$, and $\tilde{N}\equiv N_{\rm dof}/(2\pi)^{3}$ with $N_{\rm dof}$ being the number of degrees of freedom. The ${\cal R}$ and ${\cal R}_{i}$ functions appearing above are

[TABLE]

More details concerning the ${\cal H}$-functions and the manner in which they appear in the dynamical equations can be found in Refs. Alqahtani et al. (2015) and Alqahtani et al. (2017a).

To the best of our knowledge, it is not possible to analytically evaluate the ${\cal H}$-functions listed above. In practice, only one integral can be done analytically and we are left with integrals over $\phi$ and $p$. Evaluation of these 2d integrals is numerically intensive, making it infeasible to evaluate them in real-time during 3+1d simulations. One might consider interpolating them, however, they are functions of 4 variables $({\boldsymbol{\alpha}},\hat{m})$ and, in practice, one must separately interpolate these five functions and all derivatives necessary. As a consequence, one quickly runs into memory limitations, even on modern computers. A more efficient technique, which does not require a great deal of memory either, is needed. In the next subsection, we present a method for doing this.

Series expansions

Since, in practice, the $\alpha$’s do not evolve too far from $\alpha_{x}=\alpha_{y}=\alpha_{z}=1$, it makes sense to expand these integrals around such an isotropic point. After expanding around an isotropic point, the angular part of the integrals become trivial and one is left only with the $p$ integral which can easily be interpolated since it is only function of the mass. Before proceeding, we note that many of the ${\cal H}$-functions are related to each other by symmetries. As a result, we will present the method for evaluating ${\cal H}_{3}$ and ${\cal H}_{3x}$ and use symmetries to find the other ${\cal H}$-functions and their derivatives necessary.

To proceed, we expand around an arbitrary isotropic point defined by $\alpha_{i}^{2}\sim\delta_{0}$ using $\alpha_{i}^{2}=\delta_{0}+\delta_{i}\,\epsilon$ where $\delta_{0}$ is the point around which the expansion is performed, and $\epsilon$ is used to keep track of the order of the expansion. Based on this, we expand ${\cal R}$ and ${\cal R}_{x}$ as

[TABLE]

where $\delta_{i}=\alpha_{i}^{2}-\delta_{0}$.

As a result, ${\cal H}_{3}({\boldsymbol{\alpha}},\hat{m})$ can be written as

[TABLE]

where $\Omega(\boldsymbol{\delta},n)$ is the angular part of the integral which is trivial to evaluate

[TABLE]

and

[TABLE]

Similarly, ${\cal H}_{3x}({\boldsymbol{\alpha}},\hat{m})$ can be written as

[TABLE]

where $\Omega_{x}(\boldsymbol{\delta},n)$ is the angular part of the integral which is trivial

[TABLE]

The derivatives of both ${\cal H}_{3}$ and ${\cal H}_{3x}$ with respect to $\alpha$’s are straightforward since the ${\cal G}(a,\hat{m},\delta_{0})$ integral is independent of ${\boldsymbol{\alpha}}$.

There are symmetries of each ${\cal H}$-function that can be used for efficiently computing all derivatives necessary, for example, ${\cal H}_{3}$ is symmetric under the exchange of $\alpha$’s, i.e.,

[TABLE]

In a similar way,

[TABLE]

Using these identities, once one calculates one of these derivatives, the other ones can be determined using symmetry arguments. For example, once $\partial{\cal H}_{3}/\partial\alpha_{x}$ is known, the other derivatives with respect to $\alpha_{y}$ and $\alpha_{z}$ are related by the exchange symmetry

[TABLE]

Unlike ${\cal H}_{3}$, for ${\cal H}_{3x}$ we have only one identity to use, so we need two derivatives, $\partial{\cal H}_{3}/\partial\alpha_{x}$ and $\partial{\cal H}_{3}/\partial\alpha_{y}$ and then Eq. (43) can be used to find $\partial{\cal H}_{3}/\partial\alpha_{z}$

[TABLE]

We now turn to the derivative with respect to the fourth argument, $\hat{m}$. The only part in ${\cal H}_{3}$ and ${\cal H}_{3x}$ that involves $\hat{m}$ is the integral ${\cal G}(a,\hat{m},\delta_{0})$. Taking its derivative gives another integral which can be easily interpolated and used

[TABLE]

So,

[TABLE]

Using the symmetries obeyed by the ${\cal H}$ functions, one can evaluate ${\cal H}_{3y}$, ${\cal H}_{3z}$, and ${\cal H}_{3B}$ and their derivatives similarly

[TABLE]

Expansion points

Finally, we must specify which value(s) of $\delta_{0}$ to use and the order of the expansion. Since the $\alpha_{i}$’s are typically in a region $0<\alpha_{i}\lesssim 3$ during the dynamical evolution, we expand around two points corresponding to $\delta_{0}=1,4$ and interpolate between these two expansions in the intermediate region. For the interpolation, we define $r_{\rm min}=1.75$ and $r_{\rm max}=1.85$ where $r=\sqrt{\alpha_{x}^{2}+\alpha_{y}^{2}+\alpha_{z}^{2}}$ and, for ${\cal H}_{3}$, use

[TABLE]

In a similar way, ${\cal H}_{3x}$ and all derivatives necessary can be calculated. In all cases, we expand up to $12^{\rm th}$ order ($n\leq 12$) which was found to reproduce the direct numerical evaluation of all ${\cal H}$-functions very well.