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

TL;DR

This paper demonstrates that a unified viscous hydrodynamic model can successfully describe flow patterns in central p+p, p+Pb, and Pb+Pb collisions at 5.02 TeV, indicating a common hydrodynamic origin across different collision systems.

Contribution

The study introduces a generalized Monte Carlo Glauber model with constituent quarks and applies viscous hydrodynamics to all three collision types with a single parameter set.

Findings

Successfully describes elliptic and triangular flow in all systems

Supports a common hydrodynamic origin of flow in high-energy collisions

Uses a unified model for different collision systems

Abstract

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.

Equations

$$\begin{array}{r}{\vec{\mathbf{x}}}_1^{(i)}={\vec{\mathbf{x}}}^{(i)}+\sqrt{3(B-{\sigma }_g^2)}({\overrightarrow{\text{upchi}}}_1^{(i)}+{\overrightarrow{\text{upchi}}}_2^{(i)})/2,\\ {\vec{\mathbf{x}}}_2^{(i)}={\vec{\mathbf{x}}}^{(i)}+\sqrt{3(B-{\sigma }_g^2)}({\overrightarrow{\text{upchi}}}_1^{(i)}-{\overrightarrow{\text{upchi}}}_2^{(i)})/2,\\ {\vec{\mathbf{x}}}_3^{(i)}={\vec{\mathbf{x}}}^{(i)}-\sqrt{3(B-{\sigma }_g^2)}{\overrightarrow{\text{upchi}}}_1^{(i)}.\end{array}$$

$$T_{NN}^{(i,j)}=\sum _{k,l=1}^3\frac{1}{4\pi {\sigma }_g^2}{e}^{-|{\vec{\mathbf{x}}}_k^{(i)}-{\vec{\mathbf{x}}}_l^{(j)}{|}^2/(4{\sigma }_g^2)},$$

$$P_{\text{collision}}^{(i,j)}=1-e^{-{\sigma }_{gg}{T}_{NN}^{(i,j)}},$$

$$s({\vec{\mathbf{x}}}_{\perp })=\kappa \sum _{n=1}^{N_w}\sum _{j=1}^3{\gamma }_j^{(n)}\frac{1}{2\pi {\sigma }_g^2}{e}^{-|{\vec{\mathbf{x}}}_{\perp }-{\vec{\mathbf{x}}}_j^{(n)}{|}^2/(2{\sigma }_g^2)},$$

$$P_{\Gamma }(\gamma )=\frac{{\gamma }^{1/(3\theta )-1}{e}^{-\gamma /\theta }}{{\theta }^k\Gamma (1/(3\theta ))},$$

Full text

I 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 heavy-ion collisions behaves more like a strongly interacting fluid, rather than a gas of weakly-interacting quarks and gluons Adcox et al. (2005); Arsene et al. (2005); Back et al. (2005); Adams et al. (2005). Some properties of this strongly interacting QCD fluid, such as the shear viscosity coefficient and the local speed of sound, have since been constrained Romatschke and Romatschke (2007); Song and Heinz (2008a); Dusling and Teaney (2008); Bazavov et al. (2009); Schenke et al. (2011); Borsanyi et al. (2014). 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 Abelev et al. (2013); Aad et al. (2013); Adare et al. (2013, 2015); Aad et al. (2015).

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 Dumitru et al. (2015); Schenke et al. (2015); Altinoluk et al. (2015); Esposito and Gyulassy (2015); Lappi et al. (2016); Schenke et al. (2016); Hagiwara et al. (2017), or non-hydrodynamic evolution, or non-standard final-state interactions Zhou et al. (2015); Orjuela Koop et al. (2015); Bozek et al. (2015); He et al. (2016); Romatschke (2015a) 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. Romatschke (2017); Spaliński (2016) for a discussion of why non-equilibrium 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 Ortona et al. (2009); Prasad et al. (2010); Werner et al. (2011); Bozek (2012), p+Pb Bozek (2012), 3He+Au Nagle et al. (2014); Schenke and Venugopalan (2014), p+Au and d+Au collisions Romatschke (2015b) as well as the hydrodynamic description of experimental data in small systems (see e.g. Refs. Kozlov et al. (2014); Habich et al. (2016) 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. Werner et al. (2014), 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.

II Model

The present study will be based on the super-hybrid-model superSONIC Habich et al. (2015); Romatschke (2015b), which combines pre-equilibrium dynamics (based on AdS/CFT van der Schee et al. (2013)) with viscous fluid dynamic evolution (Romatschke and Romatschke (2007); Luzum and Romatschke (2008)) and late-stage hadronic rescatterings (using the hadronic cascade code B3D Novak et al. (2014)). The superSONIC model has been used in the past to successfully predict experimental flow results in p+Au, d+Au and 3He+Au collisions at $\sqrt{s}=0.2$ TeV Aidala et al. (2016). 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 much more sophisticated models for nucleon substructure have previously been discussed in the literature, cf. Refs. Avsar et al. (2011); Yan and Ollitrault (2014); Schenke and Venugopalan (2014); McGlinchey et al. (2016); Bożek et al. (2016). Given that hydrodynamics efficiently dampens small-scale fluctuations Noronha-Hostler et al. (2016), 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.

II.1 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. Welsh et al. (2016). This model shall be referred to as the OSU (Ohio State) model, and is used to generate event-by-event initial conditions as follows:

1. (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_{\textnormal{cutoff}}^{2}]$, where $b_{\textnormal{cutoff}}$ is chosen to be large enough that the probability of a collision occurring is negligible when $b>b_{\textnormal{cutoff}}$. For our simulations, we set $b_{\textnormal{cutoff}}$ to 8, 20, and 40 fm for p+p, p+Pb, and Pb+Pb collisions, respectively. 2. (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 Vries et al. (1987). The transverse positions of the nucleons are then found by projecting their 3D positions onto the transverse plane. 3. (iii)

Each nucleon in each nucleus is divided into three valence quarks. For a given nucleon $i$, the transverse positions $\vec{\mathbf{x}}_{j}^{(i)}$ ($j=1,2,3$) of the nucleon’s valence quarks are generated in the following manner. First, two vectors $\vec{\bm{\upchi}}_{1}^{(i)}$ and $\vec{\bm{\upchi}}_{2}^{(i)}$ are sampled from Gaussian distributions, for which the respective probability density functions are $P_{1}(\vec{\bm{\upchi}}_{1})=\frac{3}{2\pi}e^{-3\vec{\bm{\upchi}}_{1}^{2}/2}$, $P_{2}(\vec{\bm{\upchi}}_{2})=\frac{1}{2\pi}e^{-\vec{\bm{\upchi}}_{2}^{2}/2}$. Then, letting $\vec{\mathbf{x}}^{(i)}$ denote the transverse position of the nucleon, the quark positions are set to

[TABLE]

Here, the parameter $B$ controls the size of the nucleons, while the parameter $\sigma_{g}$ is a measure of the width of the Gaussian-shaped gluon cloud surrounding each valence quark Welsh et al. (2016). We set $B=(0.52\textnormal{ fm})^{2}$ and $\sigma_{g}=0.46\textnormal{ fm}$. 4. (iv)

For all possible pairs $(i,j)$ consisting of a beam nucleon $i$ and a target nucleon $j$, a nucleon-nucleon overlap $T_{NN}^{(i,j)}$ is computed as

[TABLE]

and the probability that a binary collision occurs between the two nucleons is determined as

[TABLE]

where $\sigma_{gg}$ is a parameter that is fixed by requiring that across a large number of generated proton+proton events, the expected probability that a collision occurs satisfies $\langle P_{\textnormal{collision}}\rangle\pi b_{\textnormal{cutoff}}^{2}=\sigma_{NN}^{\textnormal{inel}}=70\textnormal{ mb}$. Now, based on the probability $P_{\textnormal{collision}}^{(i,j)}$, 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 Mitchell et al. (2016). 5. (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(\vec{\mathbf{x}}_{\perp})$ at a point $\vec{\mathbf{x}}_{\perp}=(x,y)$ in the transverse plane of the collision is then given by

[TABLE]

where $n$ is an index over all the wounded nucleons and $N_{w}$ is the total number of wounded nucleons. $\kappa$ is a parameter whose value is chosen to reproduce final-state charged hadron multiplicities, and $\gamma^{(n)}_{j}$ is proportional to the amount of entropy deposited near midrapidity by the $j^{\textnormal{th}}$ quark of the $n^{\textnormal{th}}$ wounded nucleon. $\gamma^{(n)}_{j}$ is allowed to fluctuate from quark to quark with a probability density function

[TABLE]

where $\theta=0.75$ Welsh et al. (2016). 6. (vi)

In the last step, the continuum entropy density profile (4) is converted to an energy density $\epsilon({\bf x}_{\perp})$ using a lattice QCD equation of state Borsanyi et al. (2014) 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\propto\kappa\sum_{n=1}^{N_{w}}\sum_{j=1}^{3}\gamma^{(i)}_{j}$. Since $dS/dy$ increases with the total multiplicity of charged hadrons produced in a collision Song and Heinz (2008b), 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 superSONIC. Examples for typical transverse energy density profiles for central p+p, p+Pb and Pb+Pb collisions are shown in Fig. 1.

II.2 superSONIC

The superSONIC model converts initial energy density profiles into spectra of identified particles that can directly be compared to experimental data (see Ref. Romatschke (2015b) for a more detailed description of the model). In brief, for each initial energy-density profile $\epsilon({\bf x}_{\perp})$, a pre-equilibrium flow profile at proper time $\tau=\sqrt{t^{2}-z^{2}}$ is generated using $\vec{v}(\tau,{\bf x})=-\frac{\tau}{3.0}\vec{\nabla}\ln\epsilon({\bf x}_{\perp})$ Vredevoogd and Pratt (2009), consistent with gauge/gravity simulations of strongly coupled matter van der Schee et al. (2013), 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 $\tau=\tau_{0}=0.25$ fm using a lattice QCD equation of state Borsanyi et al. (2014) and shear and bulk viscosity values of $\frac{\eta}{s}=0.08$ and $\frac{\zeta}{s}=0.01$, respectively. Bulk viscous effects on particle spectra are at present poorly understood Monnai and Hirano (2009) so only effects of bulk viscosity on the hydrodynamic evolution is included, and for simplicity bulk and shear relaxation times are identical, $\tau_{\Pi}=\tau_{\pi}$ Kanitscheider and Skenderis (2009). The corresponding shear viscous relaxation time is varied between $\tau_{\pi}=4\frac{\eta}{sT}$ and $\tau_{\pi}=6\frac{\eta}{sT}$ in order to quantify the sensitivity of results to non-hydrodynamic modes Romatschke (2017), where $T$ denotes the local effective temperature of the system. Large variations of observables with $\tau_{\pi}$ 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\times 100$ grid points, with lattice spacings adapted to the individual size of the collision system (cf. Fig. 1). In addition, test simulations with $200\times 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 low-temperature hadronic cascade evolution with B3D Novak et al. (2014). 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 experimental measurements at mid-rapidity. The source code to superSONIC is publicly available superSONIC (version 2.1, 2016).

III Results

Using superSONIC with OSU initial conditions for the nucleon, central p+p, p+Pb and Pb+Pb collisions at $\sqrt{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 Khachatryan et al. (2017); Aad et al. (2014); Adam et al. (2016a); Aaboud et al. (2016); Chatrchyan et al. (2013). The size of the bands shown for superSONIC calculations includes statistical errors for the simulations as well as systematic uncertainties obtained from changing the second-order transport parameter $\tau_{\pi}$. 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 $\sqrt{s}=13$ TeV, more than twice the simulated collision energy of $\sqrt{s}=5.02$ TeV.

The case of p+p collision at $\sqrt{s}=5.02$ TeV has moreover been studied as a function of multiplicity, and results for the multiplicity, mean pion transverse momentum, 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 $\langle p_{T}\rangle$ values exceed the results measured by ALICE (at $\sqrt{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 $\langle p_{T}\rangle$ results, cf. Refs. Habich et al. (2016); Monnai and Hirano (2009); Ryu et al. (2015). Given the extreme sensitivity of $\langle p_{T}\rangle$ on bulk viscosity for proton+proton collisions Habich et al. (2016), it is quite possible that including bulk corrections to spectra and/or fine tuning can lead to quantitative agreement of simulation and experiment for $\langle p_{T}\rangle$ 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 non-trivial 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 $\sqrt{s}=5.02$ TeV, the model results for $\frac{dN}{dy}$ for the 0-5% highest multiplicity events are within five percent of the experimental values at midrapidity Aad et al. (2016); Adam et al. (2016b) when converting superSONIC multiplicities to pseudorapidity distributions as $\frac{dN}{dy}=1.1\frac{dN}{d\eta}$.

IV Conclusions

Relativistic p+p, p+Pb and Pb+Pb collisions at $\sqrt{s}=5.02$ TeV and small impact parameter have been simulated event-by-event using the super-hybrid-model superSONIC. 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.

Acknowledgements.

V Acknowledgements

This work was supported in part by the Department of Energy, DOE award No DE-SC0008132 and by the UROP program at the University of Colorado Boulder. We would like to thank M. Floris, M. Habich, J. Nagle and D. Perepelitsa for fruitful discussions.