Testing Hydrodynamic Descriptions of of p+p Collisions at \sqrt{s}=7 TeV

In high energy collisions of heavy-ions, experimental findings of collective flow are customarily associated with the presence of a thermalized medium expanding according to the laws of hydrodynamics. Recently, the ATLAS, CMS and ALICE experiments found signals of the same type and magnitude in ultrarelativistic proton-proton collisions. In this study, the state-of-the-art hydrodynamic model SONIC is used to simulate the systems created in p+p collisions. By varying the size of the second-order transport coefficients, the range of applicability of hydrodynamics itself to the systems created in p+p collisions is quantified. It is found that hydrodynamics can give quantitatively reliable results for the particle spectra and the elliptic momentum anisotropy coefficient $v_2$. Using a simple geometric model of the proton based on the elastic form factor leads to results of similar type and magnitude to those found in experiment when allowing for a small bulk viscosity coefficient.

In high energy collisions of heavy-ions, experimental findings of collective flow are customarily associated with the presence of a thermalized medium expanding according to the laws of hydrodynamics. Recently, the ATLAS, CMS and ALICE experiments found signals of the same type and magnitude in ultrarelativistic proton-proton collisions. In this study, the state-of-the-art hydrodynamic model SONIC is used to simulate the systems created in p+p collisions. By varying the size of the second-order transport coefficients, the range of applicability of hydrodynamics itself to the systems created in p+p collisions is quantified. It is found that hydrodynamics can give quantitatively reliable results for the particle spectra and the elliptic momentum anisotropy coefficient v2.
Using a simple geometric model of the proton based on the elastic form factor leads to results of similar type and magnitude to those found in experiment when allowing for a small bulk viscosity coefficient.
The experimental heavy-ion program at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) has provided strong evidence for the creation of an equilibrated state of matter in ultrarelativistic collisions of heavy-ions such as gold or lead [1][2][3][4][5][6][7]. Comparing the wealth of experimental data available over a large range of collision energies to theoretical model calculations, the current consensus in the field is that the matter created in ultrarelativistic heavy-ion collisions behaves like an almost ideal fluid with very low shear viscosity over entropy ratio [8][9][10][11][12][13]. This form of matter has been dubbed the 'quark-gluon plasma'.
Only a few years ago, there was a similar consensus in the field that the systems created in proton-nucleus collisions (or d+Au collisions in the case of RHIC) did not equilibrate to form a quark gluon plasma because these systems were too small, too short-lived, and contained too few particles to behave collectively. In fact, experimental data from these light-on-heavy ion collisions was regarded as a reference system in which the quark-gluon plasma component was 'known' to be absent. Similarly, the notion that quark-gluon plasmas could be formed in high energy proton-proton collisions was mostly regarded as preposterous: how could a system consisting of a handful of particles behave as a fluid?
The consensus in the field was severely challenged, if not shattered, when experimental data for anisotropic collective flow in p+Pb, p+Au, d+Au, 3 He+Au, and most recently in proton-proton collisions became available [14][15][16][17][18]. In all of these small systems, the experimental signals turned out to be similar in type and magnitude to those found in heavy-ion collisions. Furthermore, the measurements could again be well described (and in some cases predicted) by theoretical hydrodynamic model calculations [19][20][21][22][23], such as the SONIC model [24].
The experimental finding of a large elliptic flow coefficient v 2 in high energy proton-proton collisions is particularly intriguing, because a large v 2 coefficient is typically indicative of a hydrodynamic phase in the system evolution [25]. Is it at all possible for hydrodynamics to quantitatively describe the real-time evolution of system with a linear dimension of less than 1 fm and an average of five to six particles per unit rapidity? What constraints would result on QCD transport coefficients such as shear and bulk viscosity? These questions provide the motivation for performing a hydrodynamic study of high-energy proton-proton collisions.
One of the key differences of the present study with respect to most previous hydrodynamic studies of proton-proton collisions such as those in Refs. [26][27][28][29] is the inclusion of both shear and bulk viscous effects in the hydrodynamic evolution. (Note that shear viscous effects were already included in Ref. [30], which will be discussed below in more detail). Another perhaps novel aspect of the present study is that 'typical' proton-proton collisions (as opposed to high-multiplicity events such as those studied in Ref. [31]) will be discussed. Finally, the main emphasis of the present study will be a quantitative test of applicability of hydrodynamics to small systems, which has never been attempted before.

I. METHODOLOGY
In the present study, we use the hydrodynamic model SONIC [24] to simulate the matter created in proton-proton collisions. SONIC simulates the dynamics in the plane transverse to the beam axis using causal relativistic hydrodynamics in the presence of shear and bulk viscosity, followed by the hadron cascade afterburner B3D [32] in the hadronic phase for temperatures T < 0.17 GeV, while assuming boost-invariance in the longitudinal direction (see Ref. [24] for a detailed discussion of SONIC's components). It should be noted that while SONIC implements shear viscous effects when switching from hydrodynamics to the hadron cascade simulation [33], the consistent implementation of bulk viscous effects on particle spectra is currently poorly understood [34]. For this reason, bulk viscous contributions to the initial particle spectra in the hadron cascade are not included in the present description. This is different from other works in the literature (e.g. [35,36]) which use a form of the bulk viscous corrections based on a quasi-particle model [37].
SONIC is known to successfully describe experimental data for p+Pb and d+Au collisions at √ s = 5.02 TeV and √ s = 0.2 TeV collision energies, respectively, and has been used to make accurate predictions for v 2 , v 3 for 3 He+Au collisions and p+Au collisions at √ s = 200 GeV [20,23]. In order to simulate proton-proton collisions, a model for the hydrodynamic initial conditions, such as the energy density distribution in the transverse plane, is needed. These initial conditions are poorly constrained from first-principles calculations, so a basic model built on the proton form factor was used, which are described below. Besides the initial conditions, the hydrodynamic evolution in SONIC requires specification of the simulated ratios of shear viscosity and bulk viscosity to entropy density, η s , ζ s , respectively. In the following, both of these ratios were taken to be constant in temperature for simplicity. Finally, SONIC requires specification of second order transport coefficients, such as the shear and bulk relaxation times τ π , τ Π , respectively (cf. Ref. [38]). For simplicity, we have set τ π = τ Π (cf. Ref. [39]).
The value of these relaxation times controls the size of second-order gradient terms in the hydrodynamic expansion. Varying the relaxation times thus allows one to quantify the importance of second-order gradient terms in final results, and thus provides a measure of the quantitative reliability of the hydrodynamic gradient expansion. The "conventional" criterion for the applicability of hydrodynamics states that the mean free path λ needs to be much smaller than the system size L. The ratio λ L is referred to as Knudsen number, and the conventional criterion quantifies the size of first-order gradient corrections (viscous effects) to ideal hydrodynamics. In recent years there has been mounting evidence from exact solutions of far-from equilibrium quantum field theories that (second-order) hydrodynamics quantitatively applies in cases where first-order (viscous) corrections to ideal hydrodynamics are large (order unity, cf. Refs. [40][41][42][43]). Thus it may be that the "conventional" Knudsen number criterion considerably underestimates the applicability of hydrodynamics. Instead, it has been suggested that the true criterion for the applicability of hydrodynamics is set by the location of the first non-hydrodynamic singularity in the complex frequency plane [44]. In second-order hydrodynamics, the location of this pole is controlled by the value(s) of the relaxation time. Hence it is plausible that varying the relaxation time τ π allows a modern, realistic, quantitative and easily implementable test for the applicability of hydrodynamics. This is consistent with the notion of large first, but small second-order hydrodynamic corrections.
It is well known that for fixed shear-viscosity over entropy ratio, the value of τ π varies very little (only by about a factor of two) when the interaction strength in a quantum field theory is changed from zero to infinity [45,46]. With this result in mind, we choose to quantify the applicability of hydrodynamics by varying the relaxation times by 50 percent around a fiducial value of τ π = 6 η sT . If the resulting variations in the final results are large, then hydrodynamics does not apply. Conversely, if the variations turn out to be small, then this provides evidence that hydrodynamics can give a quantitatively reliable description of the system.
Basic model for the proton: We consider the initial transverse energy density distribution ε to be given by where x ⊥ = (x, y) are the coordinates in the transverse plane, τ 0 is the initialization time of hydrodynamics, b is the impact parameter of the collision, κ(τ 0 ) is an overall normalization that is fixed by the experimental multiplicity in minimum-bias collisions, and T 1,2 is the transverse charge density distribution of proton 1 and 2, respectively. The expert reader will recognize Eq. (1) as an optical-Glauber model for protons, where it should be pointed out that for protons the binary collision scaling coincides with the number of participants scaling because A = 1. Indeed, in the basic initial condition model (referred to as 'RND' for 'round' in the following), we take T (x, y) to be given by the Fourier-transform of the proton form factor F (Q 2 ), where we take the parametrization of the form factor from Ref. [47]. In the RND model, the proton is always round, and initial conditions for ε are generated by Monte-Carlo sampling of impact parameters b ∈ [0, b max ], where the upper limit b max = 1.6 fm corresponds to approximately twice the proton radius. In a variation of the 'RND' model for initial conditions, referred to as 'FLC' for 'fluctuating' in the following, spin fluctuations of the proton are considered. Using the model from Ref. [48], the overlap function is defined as where r = (x, y, z), r = |r| and N = 4π is a normalization to ensure that protons have electric charge of unity for arbitrary unit vectorsŝ,n. In the FLC model, the proton's shape may fluctuate event-by-event, and initial conditions for ε are generated by Monte-Carlo sampling of the two unit vectorsŝ,n as well as the impact parameter of the collision b ∈ [0, b max ].

II. RESULTS
Using the basic model of the proton described in the previous section, the hydrodynamic plus cascade model SONIC was initialized at τ 0 = 0.25 fm/c and results for particle spectra and momentum anisotropies were obtained that can be directly compared to experimental measurements (cf. [24]). In Fig. 1, results for the multiplicity of unidentified charged hadrons 1 and mean pion transverse momentum are shown for proton-proton collisions at √ s = 7 TeV. The   [49]) and SONIC simulations for proton models based on the proton form factor. The error bars for the SONIC simulations include systematic uncertainties for the applicability of hydrodynamics obtained from varying second-order transport coefficients; as can be seen, those error bars are significant for neither the multiplicity nor the pion < pT >, thus indicating robust applicability of hydrodynamics for these quantities. Note that the 'RND' model has been run with different shear and bulk viscosities. While the effect of changing the shear viscosity on the multiplicity and transverse momentum is minor (not shown), even a very small bulk viscosity has a large effect on the final pion transverse momentum. multiplicity in the 40-50 percent centrality class obtained by ALICE [49] was used to set the overall constant κ in the SONIC simulations. The error bars shown for the SONIC results include the systematic uncertainties for the applicability of hydrodynamics obtained from varying second-order transport coefficients, as described above. From Fig. 1 it becomes apparent that systematic uncertainties of hydrodynamics for the particle multiplicity and mean transverse momentum are small, providing evidence that a hydrodynamic description of these quantities is feasible for proton-proton collisions. The centrality dependence of multiplicity in SONIC is broadly consistent with the experimental measurements from ALICE, with a level of disagreement that can be expected given the simplicity of the initial conditions used. Considering the mean transverse pion momentum, Fig. 1 indicates that SONIC results are extremely sensitive to the presence of bulk viscosity, as is apparent from comparing the 'RND' model results for ζ s = 0 and ζ s = 0.02. This effect originates from the modification of the fluid flow from bulk viscosity, and thus is expected to be a robust feature irrespective of the hadronization prescription used (see also the discussion in Appendix A). For the proton models used, a minimum non-zero value of ζ s was needed to bring any of the theory calculation close to the experimental data from the ALICE experiment [49] for pion mean transverse momentum. Because of the crudeness of the proton model, no effort has been made to tune transport coefficient in order to match the experimental data.
Comparisons of identified particle spectra for mid-central collisions to minimum-bias experimental data are shown in Fig.2. Again, one observes reasonable overall agreement between simulations and experiment except for the case when bulk viscosity was set to zero.
The qualitative effect of bulk viscosity reducing the mean particle momenta was observed before in heavy-ion collisions, e.g. in [34,36]. However, the effect of including the bulk viscosity in proton-proton collisions is much more pronounced than in heavy-ion collisions. Specifically, we find a factor two decrease in pion momentum originating from a bulk viscosity coefficient of ζ s = 0.02, while Ref. [36] found approximately 25 percent reduction for a bulk viscosity coefficient peaking at ζ s = 0.3 (note that such high values would likely cause cavitation in the fluid [51][52][53]). In Fig. 3, the momentum anisotropy coefficient v 2 for unidentified charged particles with p T > 0.5 GeV from SONIC, including the estimated systematic uncertainty from the hydrodynamic gradient expansion is shown. (Note that v 2 is considerably smaller when a smaller p T cut is used, cf. Ref. [54]). One finds that the systematic uncertainty originating from higher order gradient terms becomes large for dN dη < ∼ 2 and η s = 0.08. Clearly, the hydrodynamic description of v 2 in this case has broken down.
However, while the systematic uncertainty originating from higher order gradient terms is sizable, it seems that hydrodynamics nevertheless is still applicable to describing v 2 in proton-proton collisions for dN dη > ∼ 2 when η s ≤ 0.08. Since this finding disagrees with an earlier prediction by one of us in Ref. [30], this point deserves further clarification. Unlike the earlier study in Ref. [30], the present study does not use hydrodynamics for temperatures below the QCD phase transition, but instead employs a hadronic cascade simulation, thus increasing overall reliability of the model. Also shown in Fig. 3 is the range of experimental results for v 2 as measured by the ATLAS experiment [18] for p+p Charged Hadron Spectra (π + +π -)/2 (K + +K -)/20 ALICE RND, η/s=0.08, ζ/s=0.00 RND, η/s=0.04, ζ/s=0.02 FLC, η/s=0.04, ζ/s=0.02 FIG. 2. Pion and kaon spectra for the 40-50 percent centrality class compared to measured minimum bias spectra for √ s = 7 TeV from the ALICE experiment [50]. The error bars for the SONIC simulations include systematic uncertainties for the applicability of hydrodynamics obtained from varying second-order transport coefficients; these error bars are smaller than the symbol size for particle spectra, thus indicating robust applicability of hydrodynamics for this quantities. Note that the 'RND' model has been run with different shear and bulk viscosities, indicating the sensitivity of particle spectra to a small bulk viscosity coefficient. TeV. The error bars for the SONIC simulations include systematic uncertainties for the applicability of hydrodynamics obtained from varying second-order transport coefficients. Right: Unintegrated momentum anisotropy for unidentified charged hadrons for the 40-50 percent centrality class compared to experimental results from ATLAS [18] and CMS [55].
collisions at √ s = 2.76 TeV and √ s = 13 TeV for N ch = 50 − 60, which roughly corresponds to the 0.5-4% centrality class (cf. [55]). The SONIC model simulation results include no or only limited event-by-event fluctuations, thereby invalidating the model results for the most central collisions ( dN dη > ∼ 10) and the most peripheral collisions ( dN dη < ∼ 1). For mid-central collisions, however, the 'RND' and 'FLC' model are broadly consistent with the magnitude of the measured v 2 coefficient by the ATLAS experiment. This finding is corroborated by the second panel in Fig. 3, where the momentum dependence of the v 2 coefficient for mid-central collisions (40-50 percent centrality class) is compared to experimental data from the ATLAS and CMS experiments.
SONIC simulation results for v 2 are sensitive to both shear and bulk viscosity coefficients, and no attempt has been made to tune the value of those coefficients in order to match the experimental data in view of the crudeness of the initial condition model. Rather, one observes that with 'typical' values for η s , ζ s the SONIC model predicts a v 2 response that is of comparable to that measured by experiment.

III. CONCLUSIONS
The hydrodynamic model SONIC was used to study proton-proton collisions at √ s = 7 TeV by employing a simple parametrization of proton based on the elastic form factor. By varying the size of the second-order transport coefficients, the applicability of hydrodynamics itself to the systems created in p+p collisions could be quantified. It was found that a hydrodynamic description of the momentum anisotropy coefficient v 2 is breaking down for dN dη < ∼ 2 when η s ≥ 0.08. Conversely, it was found that hydrodynamics can give quantitatively reliable results for the particle spectra and the elliptic momentum anisotropy coefficient v 2 when dN dη > ∼ 2. While it is somewhat surprising that hydrodynamics applies even for such low multiplicities, this finding is qualitatively in line with recent results for proton-nucleus collisions in Ref. [23]. In Ref. [23] it was found that a hydrodynamic description of v 2 was found to be reliable whereas hydrodynamics would break down sequentially starting from the higher order momentum anisotropies (first v 5 , then v 4 , etc). The finding that hydrodynamics can be applied to proton-proton collisions is also consistent with recent results from gauge/gravity duality simulations in Ref. [43]. This surprising applicability of hydrodynamics to small systems becomes somewhat less mysterious if one abandons the traditional idea of a handful of quarks and gluons forming a fluid in favor of delocalized and strongly interacting fields forming a plasma. Since hydrodynamics can be derived from a gradient expansion of quantum field theory without ever employing the concept of quasi-particles [45,56], it is perfectly reasonable to expect a tiny droplet of deconfined and strongly interacting QCD matter to behave hydrodynamically, even if this droplet will eventually hadronize into only a handful of hadrons. In principle, this notion could even offer a new interpretation of the apparently thermalized particle spectra seen to e + +e − collisions.
In the context of a hydrodynamic description, the present study provided evidence that final particle mean transverse momenta in p+p collisions are strongly sensitive to the bulk viscosity coefficient. A non-vanishing minimum value of ζ s was required to match experimental measurements of mean transverse momentum. This could indicate a possible experimental path to determining the bulk viscosity coefficient in QCD. Finally, it was found that typical elliptic momentum anisotropy coefficients v 2 obtained in the hydrodynamic model are of the same magnitude as those measured by experiment.
Clearly, many aspect of the present hydrodynamic study could and should be improved when aiming at a detailed description of experimental data in the future, such as the inclusion of more realistic event-by-event fluctuations for the proton shape, or pre-equilibrium flow. However, we do not expect these future improvements of the treatment of initial conditions to affect the applicability of hydrodynamics.
To conclude, our study provides evidence that the experimental results obtained in high energy proton-proton collisions can be understood both qualitatively and quantitatively in terms of a hydrodynamic model similar to that used in heavy-ion collisions. While the present hydrodynamic model does not describe details of the experimental measurements, it is likely that more sophisticated parametrizations of the proton could bring the same level of agreement to proton-proton collisions as is now routinely seen in heavy-ion collisions. This implies that an interpretation of the formation of a quark-gluon plasma in proton-proton collisions is consistent with the experimental data, yet does not imply that it is the only such consistent interpretation. Future work is needed to improve our qualitative and quantitative understanding of these fascinating system that link the fields of high energy and nuclear physics.  In the main text, it was mentioned that bulk viscosity affects the hydrodynamic flow pattern directly. In this appendix, the effect of bulk viscosity on the temperature and fluid velocity evolution are demonstrated through snapshots during the system evolution for an 'RND' proton collision at small impact parameter (0-10 percent centrality class), shown in Fig. 4. The panels in the figure show that adding bulk viscosity changes the hydrodynamic evolution through reducing the local fluid velocity and slowing down the temperature decrease. Since particles are sampled from the local fluid cells, smaller velocities imply smaller particle momenta, which is consistent with the finding in the main text.