One fluid to rule them all: viscous hydrodynamic description of event-by-event central p+p, p+Pb and Pb+Pb collisions at $\sqrt{s}=5.02$ TeV

The matter created in central p+p, p+Pb and Pb+Pb collisions at sqrt{s}=5.02 TeV is simulated event-by-event using the superSONIC model, which combines pre-equilibrium flow, viscous hydrodynamic evolution and late-stage hadronic rescatterings. Employing a generalization of the Monte Carlo Glauber model where each nucleon possesses three constituent quarks, superSONIC describes the experimentally measured elliptic and triangular flow at central rapidity in all systems using a single choice for the fluid parameters, such as shear and bulk viscosities. This suggests a common hydrodynamic origin of the experimentally observed flow patterns in all high energy nuclear collisions, including p+p.


INTRODUCTION
What are the properties of the matter created in ultrarelativistic ion collisions? Obtaining an answer to this question has been one of the key goals of the high energy nuclear physics community and the driving force behind the experimental heavy-ion program at both the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). Much progress has been made, such as the realization that the matter created in heavyion collisions behaves more like a strongly interacting fluid, rather than a gas of weakly-interacting quarks and gluons [1][2][3][4]. Some properties of this strongly interacting QCD fluid, such as the shear viscosity coefficient and the local speed of sound, have since been constrained [5][6][7][8][9][10]. Others, such as the minimum possible size for a QCD fluid droplet, remain yet to be unambiguously determined. Before the present decade, the mainstream expectation was that a strongly interacting QCD fluid could only be formed in "large" systems, such as those created in heavy-ion collisions. "Small" systems, such as those formed in proton+nucleus or proton+proton collisions, were not expected to flow. It thus came as a surprise to many when experimental data from proton+nucleus and proton+proton collisions both at RHIC and the LHC unambiguously demonstrated the existence of flow in these small systems [11][12][13][14][15].
The focus of the theory community has since shifted towards understanding the origin of experimentally observed flow signals in small systems. At present, there are two main schools of thought. One maintains that while experimental evidence leaves no doubt that there are flow-like signals in small systems, these signals are unrelated to those observed in heavy-ion collisions and are caused by either initial-state correlations [16][17][18][19][20][21][22], or non-hydrodynamic evolution, or non-standard final-state interactions [23][24][25][26][27] or a combination of these. The other school of thought, on which the present work will be based, adheres to Heraclit's principle of "Panta Rhei" ("Everything Flows"). According to Panta Rhei, there is no fundamental difference between the experimental flow signals in small and large systems, and both can be quantitatively explained using the laws of hydrodynamics. (See Refs. [28,29] for a discussion of why nonequilibrium hydrodynamics may be applicable to QCD fluid droplets as small as 0.15 fm). Previous work on this subject includes the prediction of flow signals in p+p [30][31][32][33], p+Pb [33], 3 He+Au [34,35], p+Au and d+Au collisions [36] as well as the hydrodynamic description of experimental data in small systems (see e.g. Refs. [37,38] for recent examples).
One of the main criticisms of the Panta Rhei approach to relativistic ion collisions has been that a hydrodynamic description of experimental data of one or two individual collision systems could be a coincidence, and that a simultaneous description of small and large systems is required. Indeed, with the notable exception of Ref. [39], previous hydrodynamic studies have focused on describing either proton+proton, proton+nucleus, or nucleus+nucleus collisions individually, rather than all of those systems together, which provides the motivation for the present study.

MODEL
The present study will be based on the superhybrid-model superSONIC [36,40], which combines preequilibrium dynamics (based on AdS/CFT [41]) with viscous fluid dynamic evolution ( [5,42]) and late-stage hadronic rescatterings (using the hadronic cascade code B3D [43]). The superSONIC model has been used in the past to successfully predict experimental flow results in p+Au, d+Au and 3 He+Au collisions at √ s = 0.2 TeV [44]. The main addition to the superSONIC model implemented here are initial conditions which allow for nucleon substructure. For the present work, a variant of the constituent quark model for the nucleon substructure will be used. This model is rather simplistic, and arXiv:1701.07145v1 [nucl-th] 25 Jan 2017 much more sophisticated models for nucleon substructure have previously been discussed in the literature, cf. Refs. [35,[45][46][47][48]. Given that hydrodynamics efficiently dampens small-scale fluctuations [49], differences between models are expected to quickly dissipate, and it is likely that other nucleon models will give almost identical results to the constituent quark model employed in this work, as long as some basic granularity of the nucleon is maintained. Hence, it is unlikely that the present work can be used to constrain nucleon structure models, with the possible exception of broad features such as the event-by-event "ellipticity" of the nucleon.

Initial Conditions
The superSONIC model requires as input initial conditions for the energy density in the transverse plane; this density should correspond to the energy deposition following the collision of two relativistic nuclei. Each nucleus consists of individual nucleons whose positions and binary collisions are treated within a Monte Carlo Glauber framework. In addition, nucleon substructure is implemented through a type of quark model of the nucleon by modifying the work of Ref. [50]. This model shall be referred to as the OSU (Ohio State) model, and is used to generate event-by-event initial conditions as follows: (i) The transverse positions of a beam and target nucleus are produced, separated by a random impact parameter of magnitude b. The impact parameter is obtained by sampling b 2 from a uniform distribution on the interval [0, b 2 cutoff ], where b cutoff is chosen to be large enough that the probability of a collision occurring is negligible when b > b cutoff . For our simulations, we set b cutoff to 8, 20, and 40 fm for p+p, p+Pb, and Pb+Pb collisions, respectively.
(ii) If either the beam or the target nucleus is simply a proton, then a single nucleon is placed at the center of that nucleus. Otherwise, if a given nucleus is a lead nucleus, then the 3D positions of the nucleons comprising the nucleus are sampled from a Woods-Saxon distribution, using a radius of 6.62 fm and a skin depth of 0.546 fm [51]. The transverse positions of the nucleons are then found by projecting their 3D positions onto the transverse plane.
(iii) Each nucleon in each nucleus is divided into three valence quarks. For a given nucleon i, the transverse positions x (i) j (j = 1, 2, 3) of the nucleon's valence quarks are generated in the following manner. First, two vectors χ (i) 1 and χ (i) 2 are sampled from Gaussian distributions, for which the respective probability density functions are P 1 ( denote the transverse position of the nucleon, the quark positions are set to (1) Here, the parameter B controls the size of the nucleons, while the parameter σ g is a measure of the width of the Gaussian-shaped gluon cloud surrounding each valence quark [50]. We set B = (0.52 fm) 2 and σ g = 0.46 fm.
(iv) For all possible pairs (i, j) consisting of a beam nucleon i and a target nucleon j, a nucleon-nucleon overlap T N N is computed as and the probability that a binary collision occurs between the two nucleons is determined as where σ gg is a parameter that is fixed by requiring that across a large number of generated pro-ton+proton events, the expected probability that a collision occurs satisfies P collision πb 2 cutoff = σ inel N N = 70 mb. Now, based on the probability P (i,j) collision , it is randomly determined whether a binary collision will occur between nucleons i and j . If a binary collision occurs, nucleons i and j are both marked as "wounded", and become participants in the collision. Only nucleons, and not individual constituent quarks are marked as "wounded", which differs from other implementations of constituent quark Monte Carlo Glauber models in the literature [52].
(v) Every wounded nucleon deposits entropy into the collision in the form of the gluon clouds surrounding its three valence quarks. (However, in the case that no nucleons are wounded, the procedure returns to step (i)). The entropy density s( x ⊥ ) at a point x ⊥ = (x, y) in the transverse plane of the collision is then given by where n is an index over all the wounded nucleons and N w is the total number of wounded nucleons. κ is a parameter whose value is chosen to reproduce final-state charged hadron multiplicities, and The actual box sizes used in simulations were adapted to individual systems. Note that for a typical p+p collision, initial conditions from the OSU model are very close to (but nevertheless slightly different from) those obtained from spherical nucleons, cf. Ref [38]. γ (n) j is proportional to the amount of entropy deposited near midrapidity by the j th quark of the n th wounded nucleon. γ (n) j is allowed to fluctuate from quark to quark with a probability density function where θ = 0.75 [50].
(vi) In the last step, the continuum entropy density profile (4) is converted to an energy density (x ⊥ ) using a lattice QCD equation of state [10] and discretized on a square lattice adapted to the size of the collision system under consideration.
Using the procedure described above, many initial energy-density profiles have been generated for p+p, p+Pb, and Pb+Pb collisions. For each initial profile, there is an associated total entropy per unit rapidity dS/dy ∝ κ Nw n=1 3 j=1 γ (i) j . Since dS/dy increases with the total multiplicity of charged hadrons produced in a collision [53], all initial density profiles are ordered into centrality classes based on their values for dS/dy. A subset of 100 initial conditions are randomly selected from each centrality class for further processing with super-SONIC. Examples for typical transverse energy density profiles for central p+p, p+Pb and Pb+Pb collisions are shown in Fig. 1. The superSONIC model converts initial energy density profiles into spectra of identified particles that can directly be compared to experimental data (see Ref. [36] for a more detailed description of the model). In brief, for each initial energy-density profile (x ⊥ ), a preequilibrium flow profile at proper time τ [54], consistent with gauge/gravity simulations of strongly coupled matter [41], while the value of the shear and bulk stress tensors will be set to zero. Using these initial conditions, 2+1 dimensional hydrodynamic simulations at mid-rapidity are then started at time τ = τ 0 = 0.25 fm using a lattice QCD equation of state [10] and shear and bulk viscosity values of η s = 0.08 and ζ s = 0.01, respectively. Bulk viscous effects on particle spectra are at present poorly understood [55] so only effects of bulk viscosity on the hydrodynamic evolution is included, and for simplicity bulk and shear relaxation times are identical, τ Π = τ π [56]. The corresponding shear viscous relaxation time is varied between τ π = 4 η sT and τ π = 6 η sT in order to quantify the sensitivity of results to non-hydrodynamic modes [28], where T denotes the local effective temperature of the system. Large variations of observables with τ π are indicative of a breakdown of hydrodynamics, while small variations suggest that hydrodynamics still applies as an effective bulk description. Simulations were performed on lattices with 100 × 100 grid points, with lattice spacings adapted to the individual size of the collision system (cf. Fig. 1). In addition, test simulations with 200 × 200 gridpoints were used to ensure that finite volume and finite resolution artifacts are under control. Once the local temperature reaches T = 0.17 GeV in a given fluid cell, hydrodynamic variables and location of the cell are stored for further processing using the lowtemperature hadronic cascade evolution with B3D [43]. B3D simulates the s-wave scatterings with a constant cross section of 10 mb and interactions through hadron resonances in the particle data book with masses up to 2.2 GeV. After resonances have stopped interacting, the final charged particle multiplicity as well as hadron spectra are obtained, and can be directly compared to exper- imental measurements at mid-rapidity. The source code to superSONIC is publicly available [57].

RESULTS
Using superSONIC with OSU initial conditions for the nucleon, central p+p, p+Pb and Pb+Pb collisions at √ s = 5.02 TeV have been simulated using one single fluid framework with fixed values of shear and bulk viscosity coefficients for all systems. The results for the differential elliptic, triangular and quadrupolar flow at midrapidity from superSONIC are shown in Fig. 2 together with experimental results from the ALICE, CMS and ATLAS experiments [58][59][60][61][62]. The size of the bands shown for su-perSONIC calculations includes statistical errors for the simulations as well as systematic uncertainties obtained from changing the second-order transport parameter τ π . The size of the uncertainty bands suggests that simulation results for all systems shown are not strongly sensitive to the presence of other, non-hydrodynamic modes, and thus a hydrodynamic effective description seems applicable.
Overall, Fig. 2 implies good agreement between the superSONIC model and experiment at low momenta for all collision systems when taking into account the systematic and statistical uncertainties in both the theory and experimental results. It should be pointed out that no fine-tuning of superSONIC parameters has been attempted, so no precision fit of the experimental data can be expected. Furthermore, note that in the case of p+p collisions, ATLAS data for v 3 , v 4 is only available for √ s = 13 TeV, more than twice the simulated collision energy of √ s = 5.02 TeV. The case of p+p collision at √ s = 5.02 TeV has moreover been studied as a function of multiplicity, and results for the multiplicity, mean pion transverse momen-tum, and integrated elliptic flow are shown in Fig. 3 together with experimental data. This figure suggests that the multiplicity distribution is well represented in the superSONIC model, while the pion mean transverse momentum only qualitatively matches experimental results: the simulated p T values exceed the results measured by ALICE (at √ s = 7 TeV) at all multiplicities. This finding is not surprising given that present simulations did not include bulk viscous corrections to the pion spectra, which can be expected to considerably affect p T results, cf. Refs. [38,55,63]. Given the extreme sensitivity of p T on bulk viscosity for proton+proton collisions [38], it is quite possible that including bulk corrections to spectra and/or fine tuning can lead to quantitative agreement of simulation and experiment for p T in p+p collisions, while not significantly altering results for p+Pb and Pb+Pb collisions. Such fine-tuning is left for future work.
Also shown in Fig. 3 is the integrated elliptic flow coefficient as a function of multiplicity, indicating that v 2 saturates at high multiplicities similar to what is observed experimentally. At low multiplicities, experimental procedures employed by different experiments lead to different results. So while the method employed by the ATLAS experiment suggests a near constant behavior of v 2 as a function of multiplicity, the method employed by CMS (not shown in Fig. 3) by construction implies that integrated v 2 decreases as multiplicity is lowered. Nevertheless, reproducing the apparent saturation of integrated v 2 at around 6 percent for high multiplicities (for which both ATLAS and CMS experiments agree on) is nontrivial for any model as this trend depends on the choice of shear viscosity and nucleon initial state parameters.
For p+Pb collisions and Pb+Pb collisions at √ s = 5.02 TeV, the model results for dN dy for the 0-5% highest multiplicity events are within five percent of the experimental values at midrapidity [64,65] when converting super- FIG. 3. Multiplicity, pion mean transverse momentum and integrated elliptic flow coefficient for p+p collisions at √ s = 5.02 TeV from superSONIC (bands) compared to experimental data from ALICE at √ s = 7 TeV and ATLAS at √ s = 5.02 TeV (symbols) [66,67]. Simulation parameters used were η s = 0.08 and ζ s = 0.01 for superSONIC and multiplicities were converted to charged particles per unit pseudorapidity as dN dy = 1.1 N ch |∆η| .

CONCLUSIONS
Relativistic p+p, p+Pb and Pb+Pb collisions at √ s = 5.02 TeV and small impact parameter have been simulated event-by-event using the super-hybrid-model su-perSONIC. Using initial conditions that allow for nucleon substructure in the form of three valence quarks, it was found that flow in all collision systems can be described simultaneously with a single set of fluid parameters. This finding suggests that the experimentally observed flow signals in proton+proton, proton+nucleus and nucleus+nucleus collisions are of common, and hydrodynamic, origin. However, more work will be needed to corroborate this conclusion.