Quasiparticle anisotropic hydrodynamics for central collisions

# Quasiparticle anisotropic hydrodynamics for central collisions

Mubarak Alqahtani    Mohammad Nopoush    Michael Strickland Department of Physics, Kent State University, Kent, OH 44242 United States
September 14, 2019
###### Abstract

We use quasiparticle anisotropic hydrodynamics to study an azimuthally-symmetric boost-invariant quark-gluon plasma including the effects of both shear and bulk viscosities. In quasiparticle anisotropic hydrodynamics, a single finite-temperature quasiparticle mass is introduced and fit to the lattice data in order to implement a realistic equation of state. We compare results obtained using the quasiparticle method with the standard method of imposing the equation of state in anisotropic hydrodynamics and viscous hydrodynamics. Using these three methods, we extract the primordial particle spectra, total number of charged particles, and average transverse momentum for various values of the shear viscosity to entropy density ratio . We find that the three methods agree well for small shear viscosity to entropy density ratio, , but differ at large . We find, in particular, that when using standard viscous hydrodynamics, the bulk-viscous correction can drive the primordial particle spectra negative at large which is clearly unphysical. Such a behavior is not seen in either anisotropic hydrodynamics approach, irrespective of the value of .

Quark-gluon plasma, Relativistic heavy-ion collisions, Anisotropic hydrodynamics, Equation of state, Quasiparticle model, Boltzmann equation
###### pacs:
12.38.Mh, 24.10.Nz, 25.75.-q, 51.10.+y, 52.27.Ny

## I Introduction

Ultrarelativistic heavy-ion collision experiments at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) study the quark-gluon plasma (QGP) using high energy nuclear collisions. The collective behavior seen in these experiments is quite successfully described by relativistic fluid dynamics. In early works, relativistic ideal hydrodynamics was applied assuming the QGP to behave like a perfect fluid Huovinen et al. (2001); Hirano and Tsuda (2002); Kolb and Heinz (2003). Later on, to include the dissipative (viscous) effects, viscous hydrodynamics has been applied Muronga (2002, 2004); Muronga and Rischke (2004); Heinz et al. (2006); Baier et al. (2006); Romatschke and Romatschke (2007); Baier et al. (2008); Dusling and Teaney (2008); Luzum and Romatschke (2008); Song and Heinz (2009); Heinz (2010); Bozek and Wyskiel (2009); Bozek (2010); El et al. (2010); Peralta-Ramos and Calzetta (2009, 2010); Denicol et al. (2010a, b); Schenke et al. (2011a, b); Bozek (2011, 2012); Niemi et al. (2011, 2012); Bożek and Wyskiel-Piekarska (2012); Denicol et al. (2012a, b); Peralta-Ramos and Calzetta (2013); Calzetta (2015); Denicol et al. (2014a); Florkowski et al. (2015a); Ryu et al. (2015); Niemi et al. (2016, 2015). Recently, due to the large momentum anisotropies generated during heavy-ion collisions, a new framework called anisotropic hydrodynamics has been developed Martinez and Strickland (2010); Florkowski and Ryblewski (2011); Ryblewski and Florkowski (2011a); Martinez and Strickland (2011); Ryblewski and Florkowski (2011b); Florkowski and Ryblewski (2012); Martinez et al. (2012); Ryblewski and Florkowski (2012); Florkowski et al. (2013a); Florkowski and Maj (2013); Ryblewski (2013); Bazow et al. (2014); Tinti and Florkowski (2014); Florkowski et al. (2014a); Florkowski and Madetko (2014); Nopoush et al. (2014); Denicol et al. (2014b); Nopoush et al. (2014); Tinti et al. (2016); Bazow et al. (2015); Nopoush et al. (2015a); Bazow et al. (2016); Alqahtani et al. (2015); Florkowski et al. (2015b); Molnar et al. (2016) (for a recent review, see Ref. Strickland (2014)). This new framework has been compared to traditional viscous hydrodynamics in many ways. For boost-invariant and transversely homogeneous systems, by comparing to exact solutions it has been shown that anisotropic hydrodynamics more accurately describes the dynamics in all cases considered Florkowski et al. (2013b, c); Bazow et al. (2014); Florkowski et al. (2014b); Denicol et al. (2014b); Tinti et al. (2016). In addition, it has been shown that anisotropic hydrodynamics best reproduces exact solutions of Boltzmann equation subject to 1+1d Gubser flow Nopoush et al. (2015b); Denicol et al. (2014c, d). Finally, we also mention that it has been shown that anisotropic hydrodynamics shows better agreement with data from ultracold Fermi gases experiments than viscous hydrodynamics Bluhm and Schaefer (2016, 2015).

The anisotropic hydrodynamics program is now focused on making phenomenological predictions for heavy-ion physics, including anisotropic freeze-out and a realistic lattice-based equation of state (EoS) Nopoush et al. (2015a). In a recent paper it was demonstrated how to impose a realistic EoS assuming approximate conformality of the QGP Nopoush et al. (2015a). In Ref. Alqahtani et al. (2015) a different method for imposing a realistic EoS was proposed in which the non-conformality of the QGP is taken into account by modeling the QGP as a gas of massive quasiparticles with temperature-dependent masses. This quasiparticle approach is motivated by perturbative results such as hard thermal loop (HTL) resummation, where the quarks and gluons can have temperature-dependent masses Weldon (1982a, b); Braaten and Pisarski (1990); Blaizot and Iancu (2002); Andersen et al. (1999, 2004); Andersen and Strickland (2005); Andersen et al. (2011); Haque et al. (2014). In the quasiparticle anisotropic hydrodynamics framework one introduces a single-finite temperature mass which is fit to available lattice data for the QCD EoS. Once is determined, the realistic EoS together with the non-equilibrium energy momentum tensor can be used to derive the dynamical equations for such a quasiparticle gas using Boltzmann equation Alqahtani et al. (2015). In this work, we extend the previous 0+1d work of Ref. Alqahtani et al. (2015) to 1+1d and we use “anisotropic Cooper-Frye freeze-out” to compute the primordial particle spectra.

Here we compare results of quasiparticle anisotropic hydrodynamics to the standard anisotropic hydrodynamics Nopoush et al. (2015a) and second-order viscous hydrodynamics. We will refer to the three methods considered herein as “aHydroQP”, “aHydro” and “vHydro”, respectively. For our comparisons, we derive the dynamical equations necessary for all cases for a 1+1d system which is azimuthally-symmetric and boost invariant. We then specialize to a relaxation-time approximation for collisional kernel and solve the partial differential equations for each approach numerically. We then present comparisons of the total number of charged particles , the average transverse momentum for pions, kaons, and protons, and the differential pion, kaon, and proton spectra predicted by each approach. We find that the three methods agree well for small shear viscosity to entropy density ratio, , but differ at large . We find, in particular, that when using standard viscous hydrodynamics, the bulk-viscous correction can drive the primordial particle spectra negative at large . Such a behavior is not seen in either anisotropic hydrodynamics approach, irrespective of the value of .

The structure of the paper is as follows. In Sec. II, we specify the notation and conventions used in the paper. In Sec. III, we review the necessary setup including basis vectors necessary in different cases, the anisotropic distribution function, and the lattice based equation of state. In Sec. IV, the Boltzmann equation and its generalization to quasiparticles with temperature-dependent masses is discussed. In Sec. V, we derive expressions for the particle four-current, energy density, and components of the pressure by taking different moments of distribution function. In Sec. VI, the dynamical equations for massive anisotropic hydrodynamics are derived for azimuthally-symmetric boost-invariant systems. In Sec. VII we discuss anisotropic Cooper-Frye freeze-out in the context of leading-order anisotropic hydrodynamics. In Sec. VIII, our numerical results obtained using the three methods for central Pb-Pb and p-Pb collisions at LHC energies are presented. Sec. IX contains our conclusions and an outlook for the future. In App. A the basis vectors used in the paper are presented. In App. B we present details about second-order viscous hydrodynamics equations. Finally, all necessary identities and function definitions are collected in Apps. C and D.

## Ii Conventions and Notation

A parentheses in the indices indicates a symmetrized form, e.g. . The metric is taken to be in the “east coast convention” such that in Minkowski space with the measure is . We also use the standard transverse projector, . When studying relativistic heavy-ion collisions, it is convenient to transform to Milne coordinates defined by , which is the longitudinal proper time, and , which is the longitudinal spacetime rapidity. For a system which is azimuthally symmetric with respect to the beam-line, it is convenient to transform to polar coordinates in the transverse plane with and . In this case, the new set of coordinates defines polar Milne coordinates. Finally, the invariant phase space integration measure is defined as

 dP≡Ndofd3p(2π)31E=~Nd3pE, (1)

where with being the number of degrees of freedom.

## Iii Setup

In this paper, we derive non-conformal anisotropic hydrodynamics equations for a system of quasiparticles with a temperature-dependent mass. To accomplish this goal, we take moments of Boltzmann equation, appropriate for the system of quasiparticles with thermal mass, to obtain the necessary anisotropic hydrodynamics equations. Using a general set of basis vectors and some simplifying assumptions relevant for 1+1d system, equations are expanded and simplified to the form appropriate for describing a boost-invariant and azimuthally-symmetric QGP. The obtained 1+1d equations are then solved numerically for our tests.

### iii.1 Ellipsoidal form

In non-conformal anisotropic hydrodynamics one introduces an anisotropy tensor of the form Nopoush et al. (2014); Martinez et al. (2012)

 Ξμν=uμuν+ξμν−ΔμνΦ, (2)

where is four-velocity, is a symmetric and traceless tensor, and is associated with the bulk degree of freedom. The quantities , , and are understood to be functions of spacetime and obey , , , and ; therefore, one has . At leading order in the anisotropic hydrodynamics expansion one assumes that the one-particle distribution function is of the form

 f(x,p)=fiso(1λ√pμΞμνpν), (3)

where has dimensions of energy and can be identified with the temperature only in the isotropic equilibrium limit ( and ).111Herein we assume that the chemical potential is zero. We note that, in practice, need not be a thermal equilibrium distribution. However, unless one expects there to be a non-thermal distribution at late times, it is appropriate to take to be a thermal equilibrium distribution function of the form , where gives Fermi-Dirac or Bose-Einstein statistics, respectively, and gives Boltzmann statistics. From here on, we assume that the distribution is of Boltzmann form, i.e. .

Since the most important viscous corrections are the diagonal components of the energy-momentum tensor, to good approximation one can assume that with and .222For a 1+1d system this is exact since one only needs the , , and components of the energy-momentum tensor. In this case, expanding the argument of the square root appearing on the right-hand side of Eq. (3) in the LRF gives

 f(x,p)=feq⎛⎜⎝1λ ⎷∑ip2iα2i+m2⎞⎟⎠, (4)

where and the scale parameters are

 αi≡(1+ξi+Φ)−1/2. (5)

Note that, for compactness, one can collect the three anisotropy parameters into vector . In the isotropic equilibrium limit, where and , one has and and, hence,

 f(x,p)=feq(ET(x)). (6)

Out of the four anisotropy and bulk parameters there are only three independent ones. In practice, we use the three variables as the dynamical anisotropy parameters since, by using Eq. (5) and the tracelessness of , one can write in terms of the anisotropy parameters, .

### iii.2 Equation of state

Herein we consider a system at finite temperature and zero chemical potential. At asymptotically high temperatures, the pressure of a gas of quarks and gluons approaches the Stefan-Boltzmann (SB) limit. At the temperatures probed in heavy-ion collisions there are important corrections to the SB limit and at low temperatures the relevant degrees of freedom change from quarks and gluons to hadrons. The standard way to determine the QGP EoS is to use non-perturbative lattice calculations. For this purpose, we use an analytic parameterization of lattice data for the QCD interaction measure (trace anomaly), , taken from the Wuppertal-Budapest collaboration Borsanyi et al. (2010). We refer the reader to the reference Alqahtani et al. (2015) for more details.

#### Method 1: Standard equation of state

In the standard approach for imposing a realistic EoS in anisotropic hydrodynamics, one derives the necessary equations in the conformal limit and exploits the conformal multiplicative factorization of the components of the energy-momentum tensor Martinez and Strickland (2010); Florkowski and Ryblewski (2011). With this method, one relies on the assumption of factorization even in the non-conformal (massive) case. Such an approach is justified by the smallness of the corrections to factorization in the massive case in the near-equilibrium limit Nopoush et al. (2015a). For details concerning this method, we refer the reader to Refs. Ryblewski and Florkowski (2012); Nopoush et al. (2015a). Although this method is relatively straightforward to implement, it is only approximate since for non-conformal systems there is no longer exact multiplicative factorization of the components of the energy-momentum tensor. This introduces a theoretical uncertainty which is difficult to quantitatively estimate.

#### Method 2: Quasiparticle equation of state

Since the standard method is only approximate, one would like to find an alternative method for imposing a realistic equation of state in an anisotropic system that can be applied for non-conformal systems. In order to accomplish this goal, we implement the realistic EoS detailed above by assuming that the QGP can be described as an ensemble of massive quasiparticles with temperature-dependent masses. As is well-known from the literature Gorenstein and Yang (1995), one cannot simply substitute temperature-dependent masses into the thermodynamic functions obtained with constant masses because this would violate thermodynamic consistency. For an equilibrium system, one can ensure thermodynamic consistency by adding a background contribution to the energy-momentum tensor, i.e.

 Tμνeq=Tμνkinetic,eq+gμνBeq, (7)

with being the additional background contribution. The kinetic contribution to the energy momentum tensor is given by

 Tμνkinetic,eq=∫dPpμpνfeq(x,p). (8)

For an equilibrium Boltzmann gas, the number and entropy densities are unchanged, while, due to the additional background contribution, the energy density and pressure are shifted by and , respectively, giving

 neq = 4π~NT3^m2eqK2(^meq), Seq = 4π~NT3^m2eq[4K2(^meq)+^meqK1(^meq)], Eeq = 4π~NT4^m2eq[3K2(^meq)+^meqK1(^meq)]+Beq, Peq = 4π~NT4^m2eqK2(^meq)−Beq, (9)

where with implicitly depending on the temperature from here on. In order to fix , one can require, for example, the thermodynamic identity

 TSeq=Eeq+Peq=T∂Peq∂T, (10)

be satisfied. Using Eqs. (9) and (10) one obtains

 dBeqdT = −12dm2dT∫dPfeq(x,p) (11) = −4π~Nm2TK1(^meq)dmdT.

If the temperature dependence of is known, then Eq. (11) can be used to determine . In practice, in order to determine , one can use the thermodynamic identity

 Eeq+Peq=TSeq=4π~NT4^m3eqK3(^meq). (12)

Using the lattice data parameterization to compute the equilibrium energy density and pressure, one can numerically solve for . We refer the reader to Ref. Alqahtani et al. (2015) for more details.

## Iv Boltzmann equation and its moments

In this paper, we derive the necessary hydrodynamical equations by taking the moments of Boltzmann equation. In this section, we introduce the Boltzmann equation and its different moments for the general case of quasiparticles with a temperature-dependent mass. Then we simplify them for the case of massless particles, when necessary, since the massless equations are used in the standard approach.

### iv.1 Boltzmann Equation

Generally, the Boltzmann equation for on-shell quasiparticles with temperature dependent mass can be written as  Jeon and Yaffe (1996); Romatschke (2012); Alqahtani et al. (2015)

 pμ∂μf+12∂im2∂i(p)f=−C[f], (13)

where is the on-shell momentum four-vector, indexes the spatial coordinates, and . In the constant mass limit, the above Boltzmann equation simplifies to

 pμ∂μf=−C[f]. (14)

The function appearing above is the collisional kernel containing all interactions involved in the dynamics. In what follows, we specialize to the case that the collisional kernel is given by the relaxation-time approximation (RTA), however, the general methods presented here can be applied to any collisional kernel. In RTA, one has

 C[f]=pμuμτeq(f−feq). (15)

In this relation, denotes the equilibrium one-particle distribution function (6) and is the relaxation time which can depend on spacetime but which we assume to be momentum-independent. To obtain a realistic model for , which is valid for massive systems, one can relate to the shear viscosity to entropy density ratio as Anderson and Witting (1974); Czyz and Florkowski (1986)

 η(T)=τeq(T)Peq(T)15κ(^meq). (16)

In this formula the function is defined as

 κ(x) ≡ x3[3x2K3(x)K2(x)−1x+K1(x)K2(x) (17) −π21−xK0(x)L−1(x)−xK1(x)L0(x)K2(x)],

where are modified Bessel functions of second kind and are modified Struve functions. Assuming that the ratio of the shear viscosity to entropy density, , is held fixed during the evolution and using the thermodynamic relation one obtains

 τeq(T)=15¯ηκ(^meq)T(1+Eeq(T)Peq(T)). (18)

Note that, , giving

 limm→0τeq(T)=5η4Peq(T). (19)

### iv.2 Moments of Boltzmann Equation

By calculating moments of Boltzmann equation one obtains evolution equations for tensors of different ranks, with the zeroth moment giving the evolution of particle four-current, the first moment giving the evolution of the energy-momentum tensor, and the second-moment describing the evolution of a particular rank three tensor. Taking the zeroth, first, and second moments of Boltzmann equation gives, respectively

 ∂μJμ = −∫dPC[f], (20) ∂μTμν = −∫dPpνC[f], (21) ∂μIμνλ−J(ν∂λ)m2 = −∫dPpνpλC[f], (22)

where the particle four-current , energy-momentum tensor , and the rank-three tensor are given by

 Jμ ≡ ∫dPpμf(x,p), (23) Tμν ≡ ∫dPpμpνf(x,p)+Bgμν, (24) Iμνλ ≡ ∫dPpμpνpλf(x,p). (25)

We note that we have introduced the non-equilibrium background field , which is the analogue of the equilibrium background in order to guarantee that the correct equilibrium limit of is obtained. In the process of the derivation one finds that, in order to write the energy momentum conservation in the form given in Eq. (21), there must be a differential equation relating and the thermal mass Alqahtani et al. (2015)

 ∂μB=−12∂μm2∫dPf(x,p). (26)

In practice, one can use (26) to write the derivative of with respect to any variable in terms of the derivative of the thermal mass times the moment of the non-equilibrium distribution function.

## V Bulk variables

In this section, the bulk variables necessary (number density, energy density, and the pressures) are calculated by taking projections of and .

### v.1 Particle four-current

The particle four-current is defined in Eq. (23). Using Eqs. (4) and (23) one has

 Jμ=(n,0)=nuμ, (27)

where and .

### v.2 Energy-momentum tensor

The energy-momentum tensor is defined in Eq. (24). Expanding it using the basis vectors one obtains

 Tμν=Euμuν+PxXμXν+PyYμYν+PzZμZν. (28)

#### v.2.1 Quasiparticle method

Using Eqs. (4), (24), and (28) and taking projections of one can obtain the energy density and the components of pressure

 E = H3(α,^m)λ4+B, Px = H3x(α,^m)λ4−B, Py = H3y(α,^m)λ4−B, Pz = H3L(α,^m)λ4−B, (29)

where .

#### v.2.2 Standard method

For a massless conformal Boltzmann gas, one has and . Using these relations, one can rewrite Eqs. (29) in the standard case, by taking the massless limit, , and hence

 E = Eeq(λ)^H3(α), Px = Peq(λ)^H3x(α), Py = Peq(λ)^H3y(α), Pz = Peq(λ)^H3L(α). (30)

These formulas suggest that, in order to impose a realistic EoS, one only has to replace and by the results obtained from lattice QCD calculations.

## Vi Dynamical equations

In order to obtain the dynamical equations from Eqs. (20)-(22), one needs to impose the RTA collisional kernel. Enforcing energy-momentum tensor conservation leads to an extra matching (constraint) equation. In total, one has the following four general dynamical equations

 ∂μJμ = 1τeq(neq−n), (31) ∂μTμν = 0, (32) ∂μIμνλ−J(ν∂λ)m2 = uμτeq(Iμνλeq−Iμνλ), (33) Ekinetic = Ekinetic,eq. (34)

Using the tensor decomposition of , , and in the basis vectors appropriate for a boost-invariant and azimuthally symmetric 1+1d system, one can expand the expressions above to obtain the final form of the equations Nopoush et al. (2014); Alqahtani et al. (2015). Choosing two equations from the first moment, three from the second moment, together with the matching condition we end up with six equations for six independent variables , , , and . The non-trivial equations from the first moment are

 DuE+Eθu+PxDxθ⊥ +1rPysinhθ⊥+1τPzcoshθ⊥=0, (35) DxPx+Pxθx+EDuθ⊥ −1rPycoshθ⊥−1τPzsinhθ⊥=0. (36)

The convective derivatives and divergences , with , are defined in App. C. Depending on the model, one can replace the bulk variables from (29) or (30) in the above equations to obtain the final form of dynamical equations for massive quasiparticle or massless standard models, respectively.

Also, taking the -, -, and - projections of the second moment equation (33), one obtains

 DuIxIx+θu+2Dxθ⊥ = 1τeq(IeqIx−1), (37) DuIyIy+θu+2rsinhθ⊥ = 1τeq(IeqIy−1), (38) DuIzIz+θu+2τcoshθ⊥ = 1τeq(IeqIz−1). (39)

Evaluating the necessary integrals using the distribution function (4), one finds

 Iu = (∑iα2i)αIeq(λ,m)+αm2neq(λ,m), (40) Ii = αα2iIeq(λ,m), (41)

with . In the massless case, they simplify to

 Iu = (∑iα2i)αIeq(λ), (42) Ii = αα2iIeq(λ), (43)

with . Finally, the matching condition is

 H3(α,^m)λ4 = H3,eq(^meq)T4, (44)

where

 H3,eq(^meq)≡limλ→Tα→1H3(α,^m). (45)

## Vii Anisotropic Freeze-out

We now turn to the topic of hadronic freeze-out. Our technique will be to perform “anisotropic Cooper-Frye freeze-out” using Eq. (3) as the form for the one-particle distribution function. This is different than the typical freeze-out prescription used in viscous hydrodynamics in which one takes into account the dissipative correction to the equilibrium distribution function only at linear order in a Taylor expansion around equilibrium. One immediate benefit of performing anisotropic freeze-out using Eq. (3) is that, with this form, one is guaranteed that the one-particle distribution function is positive-definite in all regions in phase space.

In practice, we start from the standard freeze-out integral

 N=∫Σd3ΣμJμ, (46)

where is the three-dimensional freeze-out hypersurface defining the boundary of the four-dimensional volume occupied by the fluid, is the surface normal vector, and is the particle four-current. Due to the presence of momentum-space anisotropies, one cannot simply use the momentum scale when defining the freeze-out hypersurface . Instead, one should use the energy density, from which one can obtain the effective freeze-out temperature where is obtained using our realistic EoS. After identifying , we follow the parametrization presented in Nopoush et al. (2015a).

Unlike Nopoush et al. (2015a), here we take into account the breaking of conformality. The only change required is in the distribution function itself. Parameterizing the particle momentum in the lab frame as

 pμ≡(m⊥coshy,p⊥cosφ,p⊥sinφ,m⊥sinhy), (47)

where , is the particle’s rapidity, and is the particle’s azimuthal angle. In order to set up the distribution function, having defined in Eq. (47), one can use Eqs. (2) and (3) to find in both the quasiparticle and standard cases

 pμΞμνpν = (1+Φ)[m⊥coshθ⊥cosh(y−ς)−p⊥sinhθ⊥cos(ϕ−φ)]2 (48) + ξx[m⊥sinhθ⊥cosh(y−ς)−p⊥coshθ⊥cos(ϕ−φ)]2 + ξzm2⊥sinh2(y−ς)+ξyp2⊥sin2(ϕ−φ)−Φm2.

Note that is Lorentz invariant, and by going to LRF, one can show that this is positive definite in any frame.

## Viii Results

We now turn to our numerical results. We present comparisons of results obtained using the dynamical equations of anisotropic hydrodynamics presented in Sec. VI and the second-order viscous hydrodynamics equations from Denicol et al. Denicol et al. (2012a, 2014a). For details about the vHydro equations solved herein we refer the reader to App. B.

##### Pb-Pb collisions:

For all results presented in this section we use smooth Glauber wounded-nucleon overlap to set the initial energy density. As our test case we consider Pb-Pb collisions with 2.76 GeV. The inelastic nucleon-nucleon scattering cross-section is taken to be 62 mb. For the aHydro results, we use 200 points in the radial direction with a lattice spacing of fm and temporal step size of 0.01 fm/c. For the vHydro results, we use 600 points in the radial direction with a lattice spacing of fm and temporal step size of 0.001 fm/c.333We found that the vHydro code was more sensitive to the spatial lattice spacing and required a smaller temporal step size for stability. In all cases, we use fourth-order Runge-Kutta integration for the temporal updates and fourth-order centered differences for the evaluation of all spatial derivatives.444Since the initial conditions considered herein are smooth, naive centered differences generally suffice. We take the central initial temperature to be MeV at fm/c and assume that the system is initially isotropic, i.e. for anisotropic hydrodynamics and for second-order viscous hydrodynamics. We take the freeze-out temperature to be 150 MeV in all cases shown. For the freeze-out we use 371 hadronic resonances (), with the masses, spins, etc. taken from the SHARE table of hadronic resonances Torrieri et al. (2005, 2006); Petran et al. (2014). We do not perform resonance feed down, hence all spectrum shown herein are primordial spectra.

In Figs. 1-3 we present our results for the primordial pion, kaon, and proton spectra produced for and 10, respectively. In each case, we have held the initial conditions fixed and only varied . In each of these figures the solid black line is the result from standard second-order viscous hydrodynamics (vHydro), the red short-dashed line is the result from standard anisotropic hydrodynamics (aHydro), and the blue long-dashed line is the result from quasiparticle anisotropic hydrodynamics (aHydroQP). As can be seen from Fig. 1, for , both aHydro approaches are in good agreement over the entire range shown with the largest differences occurring at low momentum. The second-order viscous hydrodynamics result, however, shows a significant downward curvature in the pion spectrum resulting in many fewer high- pions. The trend is the same for the kaon and proton spectra, however, for the larger mass hadrons the downturn is less severe. We note that although we plot only up to GeV, these plots can be extended to larger , in which case one finds that eventually the primordial pion spectrum predicted by vHydro becomes negative, which is clearly unphysical. This can be seen more clearly in Figs. 2 and 3 which show the same results for and 10, respectively. As these figures demonstrate, for larger the vHydro primordial particle spectra become unphysical at lower momenta. For example, for the differential pion spectrum goes negative at GeV while for it goes negative at GeV.

This behavior is a result of the bulk-viscous correction to the one-particle distribution function specified in Eq. (74). We have checked that if we neglect the bulk-viscous correction to the distribution function, then the resulting spectra are positive definite in the range of shown in Figs. 1-3. We have verified that this is a known issue with the bulk-viscous correction in the second-order viscous hydrodynamics approach. The same form for the bulk correction (74) was used by the authors of Ref. Ryu et al. (2015) and they also observed a downward curvature turning into negatively-valued spectra at large Denicol and Schenke (2016). This problem does not occur in either aHydro approach because the one-particle distribution function is positive-definite by construction in this framework.

Next, we turn to a discussion of Fig. 4. In this figure we plot the total number of ’s (top), ’s (middle), and ’s (bottom) obtained by integrating the differential yields over transverse momentum as a function of . The line styles are the same as in Figs. 1-3. From this figure we see that at small all approaches are in agreement, however, at large the three methods can give dramatically different results. We, in particular, note that the number of pions from vHydro drops much quicker than the two aHydro approaches. This is primarily due to the fact that in vHydro one sees a negative number of pions at large . As shown in Fig. 4 (top), the total number of pions predicted by vHydro goes negative at . This is a signal of the complete breakdown of viscous hydrodynamics and should come as no surprise since this approach is intended to be applicable only in the case of small and . Finally, we note that although aHydro and aHydroQP are in reasonably good agreement for the pion and kaon spectra, we see a rather strong dependence on the way the EoS is implemented in the proton spectra and, hence, the total number of primordial protons.

As an additional way to compare the three methods considered, in Fig. 5 we present the average for ’s (top), ’s (middle), and ’s (bottom). The line styles are the same as in Figs. 1-3. From this figure we see that both aHydro approaches predict a weak dependence of the pion and kaon on the assumed value of , whereas vHydro predicts a much more steep decrease in for the pions and kaons. Once again, we see that the vHydro for pions becomes negative for . This rapid decrease stems directly from the negativity of the pion spectra at high which is more important in this case since the integrand is more sensitive to the high- part of the spectra. Finally, we note that once again the proton spectra and hence for protons is sensitive to the way in which the EoS is implemented when comparing the two aHydro approaches.

Finally, in Fig. 6 we plot the total number of charged particles as a function of predicted by each of the three approaches considered. Once again the line styles are the same as in Figs. 1-3. As this figure demonstrates, for small all three frameworks are in agreement and approach the ideal result as tends to zero. All three frameworks predict that at first increases, then reaches a maximum, and then begins to decrease. The precise turnover point depends on the method with the lowest turnover seen using vHydro around ; however, this turnover is due in large part to the fact that the total pion number drops precipitously in vHydro, eventually becoming negative, as we pointed out in the discussion of Fig. 4. As a result, one cannot trust the large predictions of vHydro.

##### p-Pb collisions:

We now turn to the case of an asymmetric collision between a proton and a nucleus. For this purpose, we use the same parameters as used for the Pb-Pb collisions considered previously, except for p-Pb we use a lower initial temperature of MeV . In addition, for the aHydro results, we use a lattice spacing of fm and for vHydro results, we use a lattice spacing of fm. The smaller lattice spacings simply reflect the smaller system size of the QGP created in a p-Pb collision.

In Figs. 7-8 we present our results for the primordial pion, kaon, and proton spectra produced for and 3, respectively in p-Pb collisions. As we can see from Fig. 7, both aHydro approaches are in a good agreement at high , however, there are some quantitative differences at low . Comparing the low difference with that seen in Pb-Pb collisions, we find that there is a larger variation in the p-Pb spectra comparing aHydro and aHydroQP. This variation is a bit worrisome since it indicates a kind of theoretical uncertainty in the aHydro approach. Importantly, however, we mention that the two aHydro results are quite different than the vHydro result. For p-Pb collisions, the vHydro result shows the same behavior seen in the Pb-Pb collisions, namely that the particle spectra goes negative at high . In Fig. 8, one can clearly see the unphysical behavior in the vHydro primordial spectra. As a result, there is even larger theoretical uncertainty associated with applications of vHydro to p-Pb collisions.

## Ix Conclusions and Outlook

In this paper, we compared three different viscous hydrodynamics approaches: aHydro, aHydroQP, and vHydro. For all three cases we included both shear and bulk-viscous effects using the relaxation time approximation scattering kernel. In the standard aHydro approach one uses the standard method for imposing an equation of state in anisotropic hydrodynamics, which is to obtain conformal equations and then break the conformality only when introducing the EoS itself. In aHydroQP, one takes into account the breaking of conformality at the outset by modeling the QGP as a quasiparticle gas with a single temperature-dependent mass which is fit to available lattice data for the EoS. Finally, for our comparisons with viscous hydrodynamics we used the formalism of Denicol et al Denicol et al. (2012a, 2014a) specialized to the case of relaxation time approximation.

For each method, we specialized to the case of 1+1d boost-invariant and azimuthally-symmetric collisions using smooth Glauber initial conditions. We specialized to this case because of the computational intensity of the aHydroQP approach which requires real-time evaluate of complicated multi-dimensional integrals which are functions of all three anisotropy parameters and the local temperature-dependent mass. Using the resulting numerical evolution, we then extracted fixed energy density freeze-out hypersurfaces in each case and implemented the scheme-appropriate freeze-out to hadrons allowing us to have an apples-to-apples comparison between the three different approaches. We found that the primordial particle spectra, total number of charged particles, and average transverse momentum predicted by the three methods agree well for small shear viscosity to entropy density ratio, , but differ at large . Our most important finding was that when using standard viscous hydrodynamics, the bulk-viscous correction can drive the primordial particle spectra negative at large . Such a behavior is not seen in either anisotropic hydrodynamics approach, irrespective of the value of .

Looking to the future, it is feasible to extend the aHydroQP approach to (3+1)d, however, this will be numerically intensive and require parallelization to implement fully. One possibility is to use polynomial fits to parametrize the various massive -functions necessary instead of evaluating them on-the-fly in the code. If this is possible, then aHydroQP could become a viable alternative to the standard method of implementing the EoS in aHydro and, since it takes into account the non-conformality of the system from the beginning, this could give us an idea of the theoretical uncertainty associated with the EoS method used in phenomenological applications. For now, this work will serve as a reference point for possible differences between the two approaches to imposing the aHydro EoS.

###### Acknowledgements.
We thank G. Denicol and B. Schenke for useful discussions. M. Strickland and M. Nopoush were supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Awards No. DE-SC0013470. M. Alqahtani was supported by a PhD fellowship from the University of Dammam.

## Appendix A Basis Vectors

One can define the general basis vectors in the lab frame (LF) by performing the Lorentz transformation necessary to go from local rest frame to the LF. The transformation required can be constructed using a longitudinal boost along the beam axis, followed by a rotation around the beam axis, and finally a transverse boost by along the -axis Florkowski and Ryblewski (2012); Martinez et al. (2012). This parametrization gives

 uμ ≡ (coshθ⊥coshϑ,sinhθ⊥cosφ,sinhθ⊥sinφ,coshθ⊥sinhϑ), Xμ ≡ (sinhθ⊥coshϑ,coshθ⊥cosφ,coshθ⊥sinφ,sinhθ⊥sinhϑ),
 Yμ ≡ (0,−sinφ,cosφ,0), Zμ ≡ (sinhϑ,0,0,coshϑ), (49)

where the three fields , , and are functions of Cartesian Milne coordinates . Introducing another parametrization by using the temporal and transverse components of flow velocity, one has

 u0=coshθ⊥,ux=u⊥cosφ,uy=u⊥sinφ,⟹uμ≡(u0coshϑ,ux,uy,u0sinhϑ),Xμ≡(u⊥coshϑ,u0uxu⊥,u0uyu⊥,u⊥sinhϑ),Yμ≡(0,−uyu⊥,uxu⊥,0),Zμ≡(sinhϑ,0,0,coshϑ), (50)

where . For a boost-invariant and azimuthally-symmetric system, one can simplify the basis vectors by identifying and where and are the spacetime rapidity and the azimuthal angle, respectively. In this case, the basis vectors (49) simplify to

 uμ = (coshθ⊥coshς,sinhθ⊥cosϕ,sinhθ⊥sinϕ,coshθ⊥sinhς), Xμ = (sinhθ⊥coshς,coshθ⊥cosϕ,coshθ⊥sinϕ,sinhθ⊥sinhς), Yμ = (0,−sinϕ,cosϕ,0), Zμ = (sinhς,0,0,coshς). (51)

## Appendix B Second-order viscous hydrodynamics

Similar to anisotropic hydrodynamics, the viscous hydrodynamics dynamical equations can be obtained by taking moments of the Boltzmann equation. The energy-momentum tensor in the case when the bulk correction is not ignored is

 Tμν=Euμuν−Δμν(P+Π)+πμν. (52)

Taking the first and second moments of Boltzmann equation one obtains Denicol et al. (2014a, 2012a)

 (E+P+Π)Duuμ = ∇μ(P+Π)−Δμν∇σπνσ+πμνDuuν, (53) DuE = −(E+P+Π)θu+πμνσμν, (54) τΠDuΠ+Π = −ζθu−δΠΠΠθu+φ1Π2+λΠππμνσμν+φ3πμνπμν, (55) τπΔμναβDuπαβ+πμν = 2ησμν+2τππ⟨μαων⟩α−δπππμνθu+φ7π⟨μαπν⟩α−τπππ⟨μασν⟩α+λπΠΠσμν+φ6Ππμν, (56)

where and are the equilibrium (isotropic) energy density and pressure, respectively, and are the shear and bulk relaxation time, respectively, and is the shear-shear-coupling transport coefficient. The various notations used are

 (57)

The non-vanishing Christoffel symbols for polar Milne coordinates are , , , and . Also, we note that for the smooth initial conditions considered herein in 1+1d, the vorticity tensor vanishes. As shown in Ref. Molnár et al. (2014), the terms , , , and appear only because the collision term is nonlinear in the single-particle distribution function. In the case of the RTA, the collision term is assumed to be linear in the single-particle distribution function and one has . In this case, the shear and bulk relaxation times, and , respectively, are equal to the microscopic relaxation time , i.e., Denicol et al. (2014a). The coefficients appearing in the equation for the bulk and shear corrections are Denicol et al. (2014a)

 ζτΠ=(13−c2s)(E+P)−29(E−3P),δΠΠτΠ=1−c2s,λΠπτΠ=13−c2s,ητπ=45P+115(E−3P),δππτπ=43,τππτπ=107,λπΠτπ=65, (58)

with being the speed of sound squared and .

### 1+1d viscous hydrodynamics equations of motion

In the boost-invariant and azimuthally-symmetric case, one has and, as a result, . In addition, for this case, the shear tensor has the following form

 πμν=⎛⎜ ⎜ ⎜⎝πττπτr00πτrπrr0000πϕϕ0000πςς⎞⎟ ⎟ ⎟⎠. (59)

In this case, expanding Eqs. (53), (54), and (55) in polar Milne coordinates one obtains six equations where only five of them are independent

 (E+P+Π)Duuτ = −(ur)2[∂τ(P+Π)−dνπντ]−uτur[∂r(P+Π)−dνπνr], (60) (E+P+Π)Duur = −uτur[∂τ(P+Π)−dνπντ]−(uτ)2[∂r(P+Π)−dνπνr], (61) DuE = −(E+P+Π)θu−πrr(1−v2)2∇⟨rur⟩−