Fluid dynamics with saturated minijet initial conditions in ultrarelativistic heavy-ion collisions

Using next-to-leading order perturbative QCD and a conjecture of saturation to suppress the production of low-energy partons, we calculate the initial energy densities and formation times for the dissipative fluid dynamical evolution of the quark-gluon plasma produced in ultrarelativistic heavy-ion collisions. We identify the framework uncertainties and demonstrate the predictive power of the approach by a good global agreement with the measured centrality dependence of charged particle multiplicities, transverse momentum spectra and elliptic flow simultaneously for the Pb+Pb collisions at the LHC and Au+Au at RHIC. In particular, the shear viscosity in the different phases of QCD matter is constrained in this new framework simultaneously by all these data.

The main goal of ultrarelativistic heavy-ion collisions at the Large Hadron Collider (LHC) and the Relativistic Heavy-Ion Collider (RHIC) is to determine the thermodynamic and kinetic properties of strongly interacting matter. The measured hadronic transverse momentum (p T ) spectra at the LHC and RHIC provide convincing evidence for a formation of a strongly collective system and a nearly thermalized quark-gluon plasma (QGP) [1]. In particular, the observed systematics of the Fourier harmonics v n = cos(nφ) of the azimuth-angle distributions, are remarkably consistent with a low-viscosity QCD matter whose expansion and cooling are describable with dissipative relativistic fluid dynamics [2][3][4][5][6][7][8][9][10][11][12].
The essential inputs to the fluid dynamics are the initial energy density and flow of the matter created in the collision. However, the final state observables like multiplicities, p T spectra and v n , are also strongly affected through the fluid dynamical expansion by the viscosity and the equation of state (EoS). Thus the entire spacetime evolution, including partons in the colliding nuclei, the primary production and thermalization of QCD matter and the subsequent fluid dynamical evolution, becomes highly convoluted. Description of all these dynamics in a coherent way, leading to quantitative predictions and a meaningful determination of the QCD matter properties from the measurements, provides an ultimate challenge in the field. As discussed in this paper, the determination of e.g. the temperature dependence of the shear viscosity-to-entropy ratio η/s(T ) calls for a simultaneous theory analysis of all possible bulk (low-p T ) observables at the LHC and RHIC.
Parton saturation is a viable mechanism to control the otherwise unsuppressed production of soft small-p T quanta in hadronic and nuclear collisions [13][14][15][16]. In essence saturation means that there exists a semihard scale controlling the particle production in the collision. In the perturbative QCD (pQCD) + saturation framework we consider here, the primary particle production in A+A collisions is computed in collinear factorization by approaching the saturation at semi-hard scales from the perturbatively controllable high-p T side [17,18]. Perturbative QCD provides an excellent description of hard processes in hadronic and nuclear collisions at interaction scales Q 1 GeV [19]. Moreover, this framework allows for a quantification of the particle production uncertainties, and their propagation through the fluid dynamical evolution in nuclear collisions [18]. In addition to the internal consistency of the pQCD-based approach, it should be noted that perturbative primary gluon production in heavy-ion collisions is complementary to the Color-Glass Condensate models [20] which build on soft gluon fields. If these different high-energy QCD approaches produce similarly successful heavy-ion phenomenology, the overall uncertainty in determining the QCD matter properties can be dramatically reduced.
The present work has roots in the so-called EKRT saturation model [17], which successfully predicted the multiplicities and p T spectra in central A+A collisions at RHIC and LHC [21][22][23][24], and also the centrality dependence at RHIC [25] (cf. Fig. 23(a) in [26]). Here we use the next-to-leading-order (NLO)-improved pQCD + saturation framework of [18] to calculate the initial QGP energy density profiles and formation times, and combine these with viscous fluid dynamics. We analyse the centrality dependence of charged particle multiplicities, p T spectra and elliptic flow (v 2 ) at the LHC and RHIC in terms of the few physical key-parameters of the framework. We show that a good simultaneous description of all these observables can indeed be obtained without retuning the framework from one collision system (cmsenergy, nuclei, centrality) to another. This results in the robust predictive power of the approach, originating from the pQCD calculation of the QGP initial conditions. Most importantly, this predictive power enables us to study and restrict the ratio η/s(T ) in the different QCD-matter phases more consistently in a simultaneous multiobservable analysis of the LHC and RHIC data.
Let us then discuss the details of our framework [18]. The rigorously calculable part is the minijet E T production in an A+A collision, in a rapidity interval ∆y and above a p T scale p 0 , where s = (x, y) is the transverse location, b the impact parameter, and T A (s) the standard nuclear thickness function with the Woods-Saxon nuclear density profile. The first E T -moment of the minijet E T distribution, σ E T p0,∆y,β , [18,27] is in NLO where dσ 2→n are the collinearly factorized minijet production cross sections and [DPS] n denote the phase-space differentials for the 2 → 2 and 2 → 3 cases [18,28]. We apply the CTEQ6M parton distribution functions (PDFs) [29] with the EPS09s impact-parameter dependent nuclear PDFs [30]. The measurement functionsS 2 andS 3 define the hard scattering in terms of the minijet transverse momenta p T,i and the cut-off scale p 0 , as well as the total minijet E T produced in ∆y: where E T,n = n i=1 Θ(y i ∈ ∆y)p T,i and Θ is the step function. These functions, analogous to the jet definitions [31], are constructed so that σ E T p0,∆y,β is a welldefined, infrared-and collinear-safe, quantity to compute. The hardness-parameter β defines the minimum E T in the interval ∆y. As discussed in [18], any β ∈ [0, 1] is acceptable for the rigorous NLO computation.
Following the new angle in formulating the minijet saturation [18], the E T production is expected to cease when the 3 → 2 and higher-order partonic processes start to dominate over the conventional 2 → 2 processes. For a central collision of identical nuclei of radii R A this leads to a transversally averaged saturation criterion with an unknown, α s -independent, proportionality constant K sat ∼ 1. Generalizing to non-zero impact parameters and localizing in the transverse coordinate plane gives where the l.h.s. is the p 0 -dependent NLO pQCD calculation defined in Eq. (1). For given K sat and β, we solve the above equation for , and obtain the total dE T /d 2 s in a mid-rapidity unit ∆y = 1 at saturation from the r.h.s. as K sat p 3 sat /π. Once the solution p sat is known, the local energy density is obtained [21,23] as where the local formation time is τ s = 1/p sat . Fig. 1 shows examples of p sat ( √ s N N , A, s, b; K sat , β) as a function of T A T A , calculated for fixed values of K sat , β and with b = 0 and three other fixed impact parameters corresponding to the centrality classes 0-5%, 20-30% and 40-50% in √ s N N = 2.76 TeV Pb+Pb collisions at the LHC and 200 GeV Au+Au at RHIC. To a very good approximation, the b and s dependence of p sat comes only through T A T A . This is due to the weak s dependence of the nPDFs near the centres of the nuclei [30]. The approximate power-law scaling behaviour seen at large T A T A can then be understood as expained in [32]. We identify two main uncertainties in mapping the pQCD + saturation calculation to an initial state for fluid dynamics: (i) The energy density given by Eq. (5) is at a time τ s = 1/p sat , i.e. different at each transverse point s, while for fluid dynamics we need the initial condition at a fixed time τ 0 . (ii) We cannot trust the pQCD calculation down to p sat → 0, but we need to set a minimum scale p min sat ≫ Λ QCD . Wherever p sat ≥ p min sat we can use the pQCD calculation, but the other regions, i.e. low density edges, need to be treated separately.
We fix a minimum saturation scale as p min sat = 1 GeV. Correspondingly, the maximum formation time in our framework is τ 0 = 1/p min sat . Then, we evaluate the energy densities from τ s (s) to τ 0 using either the Bjorken free streaming ε(τ 0 ) = ε(τ s )(τ s /τ 0 ) (FS) or the Bjorken hydrodynamic scaling solution ε(τ 0 ) = ε(τ s )(τ s /τ 0 ) (4/3) (BJ). We take these two limits to represent the uncertainty in the early pre-thermalization evolution: In the free streaming case the transverse energy is preserved, while the other limit corresponds to the case where a maximum amount of the transverse energy is reduced by the longitudinal pressure.
To obtain the energy density ε(s, τ 0 ) in the transverse region where p sat < p min sat , we use an interpolation ε = C(T A T A ) n , where the power n = 1 2 [(k + 1) + (k − 1) tanh({σ N N T A T A − g}/δ)] with the total inelastic nucleon-nucleon cross-section σ N N , and g = δ = 0.5 fm −2 . This smoothly connects the FS/BJevolved pQCD energy density ε(p min sat ) = C(T A T A ) k to the binary profile ε ∝ T A T A at the dilute edge.
For the fluid-dynamical evolution, we use the state-ofthe art 2+1 D setup previously employed in Ref. [11,12,33], assuming longitudinal boost invariance, a zero netbaryon density and thermalization at τ 0 . The equations of motion are given by the conservation laws for energy and momentum, ∂ µ T µν = 0. The evolution equation of the shear-stress tensor π µν = T µν is given by transient relativistic fluid dynamics [34][35][36], where the co-moving time derivative u µ ∂ µ is denoted by the dot, η is the shear viscosity coefficient, σ µν = ∂ µ u ν is the shear tensor, θ = ∂ µ u µ is the expansion rate, and the angular brackets denote the symmetrized and traceless projection, orthogonal to the fluid four-velocity u µ . The coefficients of the non-linear terms are taken to be c 1 = 4τ π /3, c 2 = 10τ π /7 and c 3 = 9/(70p), where p is the thermodynamic pressure and τ π = 5η/(ε + p). For details of the numerical algorithm, see Refs. [11,37]. The hadron spectra are calculated with the Cooper-Frye freeze-out procedure [38] by using Israel's and Stewart's 14-moment ansatz for the dissipative correction to the local equilibrium distribution function, ing different hadron species and p µ i the 4-momentum of the corresponding hadron. The freeze-out temperature is here always T dec = 100 MeV. After calculating the thermal spectra, we include the contribution from all 2-and 3-particle decays of unstable resonances in the EoS.
We use the lattice QCD and hadron resonance gas (HRG) based EoS s95p-PCE-v1 [39] with a chemical freeze-out temperature T chem = 175 MeV. Although the rather high T chem leads to an overabundance of protons, it however reproduces the low-p T region of the p T -spectra much better than e.g. T chem = 150 MeV.
For a rough but realistic (non-constant [40]) shear viscosity description, we assume the ratio η/s to decrease linearly as a function of temperature in the hadronic phase, be in a minimum at the matching-temperature 180 MeV of the HRG/QGP phases in the used EoS, and either to increase or stay constant vs. T in the QGP phase [11,12]. Fig 2 shows the η/s(T ) which in our framework best reproduce the v 2 coefficients simultaneously at RHIC and LHC.
At this point, we have a fixed framework with four correlated unknowns, {K sat , β, BJ/FS, η/s(T )}, to be de-termined using the LHC and RHIC data on the centrality dependence of the charged particle multiplicities, p T spectra and v 2 . We proceed by scanning the parameters K sat = O(1), β ∈ [0, 1] and η/s(T ). In particular, we vary the minimum value and slopes of η/s(T ), keeping its general shape as in Fig. 2. Both the BJ and FS prethermal evolutions are considered. In practice, for each fixed {β, BJ/FS, η/s(T )}, the remaining parameter K sat is always tuned such that the multiplicity in the 0 − 5 % most central collisions at the LHC is reproduced. In Fig. 3a we show the computed centrality dependence of the charged hadron multiplicity in Pb+Pb collisions at √ s N N = 2.76 TeV compared with the AL-ICE data [41]. As demonstrated here, several sets {K sat , β, BJ/FS, η/s(T )} give a good agreement with the measurement. However, the data clearly favours β ∼ 1 and slightly the FS scenario over the BJ. For comparison, we also show the results obtained with the usual (non-saturation) eBC and eWN Glauber model initial states [42].
In Fig. 3b we show the multiplicities for Au+Au collisions at √ s N N = 200 GeV, using the same parameter sets {K sat , β, BJ/FS, η/s(T )} as in panel 3a, and compare with the PHENIX [43] and STAR [26] data. We note that although the RHIC data would seem to favour a slightly smaller β and the BJ case, the overall simultaneous agreement at RHIC and LHC is rather good.
As long as the centrality dependence of the multiplicity is described, all the scenarios studied here give a very good description of the charged hadron p T -spectra. More relevant parameters in this case are T chem and T dec which here are kept unchanged from RHIC to LHC. The obtained p T spectra are shown in Fig. 3c for the LHC and in Fig. 3d for RHIC. The data are from Refs. [44] and Ref. [45,46], correspondingly.
In Figs. 3e and 3f we show the elliptic flow coefficients v 2 (p T ) at the LHC and RHIC, respectively. The data are from ALICE [47] and STAR [48]. The v 2 (p T ) coefficients depend strongly on the η/s parametrization, and   e.g. an ideal fluid description (not shown) does not give a correct v 2 (p T ). By scanning the η/s(T ) as explained above, while keeping K sat of order 1, we observed that a good simultaneous agreement with the measurements is obtained with the cases shown in Fig. 2. We emphasize that at RHIC, where the flow gradients are larger, one should require the agreement in particular in the smallp T region, where the dissipative corrections to the particle distributions do not grow unphysically large. Note especially that since η/s(T ) is considered as a material property, it must not be changed between different collision systems.
To conclude, we computed the energy density profiles and formation times of the produced QGP at the LHC and RHIC in a new NLO-improved pQCD + local saturation framework of considerable predictive power. The subsequent evolution of these initial conditions was described with dissipative fluid dynamics. Identifying the framework uncertainties, a good global agreement with the measured centrality dependence of the low-p T bulk observables was obtained simultaneously at the LHC and RHIC. In particular, we were able to constrain the η/s(T ) parametrization simultaneously by all these data. In the future, we will extend this analysis to include event-byevent fluctuations.